DYNAMICAL  EVOLUTION  OF  ASTEROIDAL  AND  COMETARY 
PARTICLES  AND  THEIR  CONTRIBUTION  TO  THE  ZODIACAL  CLOUD 


BY 
JER-CHYI  LIOU 


A  DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 

OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 

OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 

DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 
1993 


jim^^'j>(:u^m>  nmm%^  m^^mo 


To  my  parents,  brother,  sister  in  law,  and  my  beloved  wife  for  their  love  and  support. 


ACKNOWLEDGMENTS 


I  would  like  to  thank  my  adviser,  Dr.  Stan  Dermott,  and  the  other  members  in 
my  committee.  Dr.  H.  Campins,  Dr.  F.  E.  Dunnam,  and  Dr.  H.  Smith  for  their 
support  and  guidance  through  out  the  work  of  this  dissertation.  Special  thanks  go  to 
Dr.  A.  Smith  for  reading  the  thesis  and  providing  me  with  useful  suggestions. 

Special  thanks  to  Dr.  Yulin  Xu  for  helping  me  understand  SIMUL  and  providing 
me  with  the  observational  data  from  IRAS.  I  would  like  to  thank  Dr.  B.  Gustafson 
for  many  useful  discussions.  I  would  also  like  to  thank  Dan  Durda  for  his  comment 
on  the  collisional  evolution  of  the  asteroids. 

Finally,  I  would  like  to  thank  my  fellow  graduate  students,  Jaydeep  Mukherjee 
for  his  help  over  the  years,  Sandra  Clements  for  showing  me  how  to  be  a  good 
observer  at  RMH  Observatory,  David  Kaufmann  and  Damo  Nair  for  making  me  the 
best  backgammon  player  in  the  department. 


m 


TABLE  OF  CONTENTS 

ACKNOWLEDGMENTS iii 

LIST  OF  TABLES vi 

LIST  OF  FIGURES vii 

ABSTRACT xv 

CHAPTERS 

1  INTRODUCTION     1 

The  Zodiacal  Cloud 1 

Objective  of  this  Dissertation 3 

2  SECULAR  PERTURBATION  THEORY    6 

Low-order  Secular  Perturbation  Theory 6 

The  Solar  System 6 

Small  Bodies  in  the  Solar  System 10 

Numerical  Approach 12 

Comparison  between  the  Analytical  and  Numerical  Calculations  .  .    14 

Sun-Jupiter-Satum  System    14 

Asteroidal  Particles  in  the  Solar  System 15 

3  ORBITAL  EVOLUTION  OF  ASTEROIDAL  PARTICLES 25 

Radiation  Pressure,  Poynting-Robertson  Drag,  and  the  Corpuscular 
Drag 25 

Static  Theory  and  Dynamical  Theory 28 

Numerical  Comparisons 32 

Drag  Effect 32 

iv 


Passage  Through  Resonance 34 

Orbital  Evolution  of  Dust  Particles  from  the  Major  Hirayama 

Asteroidal  Families 36 

9  fim  Diameter  Particles 38 

4  fim  and  14  ^m  Diameter  Particles 41 

Near-Earth  Dust  Grains 42 

4  ORBITAL  EVOLUTION  OF  COMETARY  PARTICLES 77 

Orbital  Evolution  of  Encke-type  Dust  Particles 77 

Gauss'  Equations 80 

Can  Comets  be  the  Source  of  the  Solar  System  Dust  Bands?  ....  81 

Distribution  of  Encke-type  Dust  Particles  in  the  Solar  System    ...  83 

5  STRUCTURE  OF  THE  ZODIACAL  CLOUD 106 

Contribution  from  Asteroidal  Dust  Particles 108 

Contribution  from  Encke-type  Dust  Particles 110 

Discussion Ill 

BIBLIOGRAPHY 129 

BIOGRAPHICAL  SKETCH 132 


LIST  OF  TABLES 


2.1:  Comparison  between  numerical  and  analytical  results  for  the  real 

Jupiter  and  Saturn 17 

2.2:  Comparison  between  numerical  and  analytical  results  for  small  m, 

e,  I,  Jupiter  and  Saturn 18 


VI 


LIST  OF  FIGURES 


2.1:  Variation  in  h,  k,  p,  and  q  of  Jupiter  (subscript  J)  with  time  obtained 

from  numerical  integration  of  the  Sun- Jupiter-Saturn  system    ...    19 

2.2:  Variation  in  h,  k,  p,  and  q  of  Saturn  (subscript  S)  with  time  obtained 

from  numerical  integration  of  the  Sun-Jupiter-Satum  system    ...    20 

2.3:               Distribution  in  inclination  space  of  Koronis-like  particles  at  t=0  and 
t=280,000  years 21 

2.4:  Distribution  in  eccentricity  space  of  Koronis-like  particles  at  t=0 

and  t=280,000  years,  where  u;  is  the  longitude  of  pericenter  ....    22 

2.5:               Variation  in  the  forced  p  and  q  over  a  period  of  300,000  years.  The 
dots  are  from  direct  numerical  integration  while  the  solid  line  is 
from  secular  perturbation  theory 23 

2.6:               Variation  in  the  forced  h  and  k  over  a  period  of  300,000  years.  The 
dots  are  from  direct  numerical  integration  while  the  solid  line  is 
from  secular  perturbation  theory 24 

3.1:               Variation  in  the  total  (subscript  t),  forced  (subscript  f),  and  proper 
(subscript  p)  eccentricities  and  longitudes  of  pericenter  from  t=0  to 
t=6t 47 

3.2:  Variation  in  the  osculating  inclination  vectors  of  9  ;/m  diameter  Eos- 

type  dust  particles  as  they  spiral  towards  the  Sun  from  1.8  AU.  Each 
panel  shows  the  vectors  at  various  mean  distances  from  the  Sun  .    48 

3.3:  Variation  in  the  osculating  eccentricity  vectors  of  9  //m  diameter  Eos- 

type  dust  particles  as  they  spiral  towards  the  Sun  from  1.8  AU.  Each 
panel  shows  the  vectors  at  various  mean  distances  from  the  Sun  .    49 


vu 


3.4:  Proper  inclination  as  a  function  of  semimajor  axis  of  9  nm 

diameter  Eos-type  dust  particles  as  they  spiral  towards  the  Sun 
from  1.8  AU.  Error  bars  are  the  dispersions  in  the  radii  from  the 
mean  radius  obtained  from  the  least  squares  fit 50 

3.5:  Comparison  of  the  forced  inclination  obtained  from  the  dynamical 

theory  (bold  line)  and  from  numerical  integration  (points)  of  9  /im 
diameter  Eos-type  dust  particles  as  they  spiral  towards  the  Sun  from 
1.8  AU.  Error  bars  are  the  dispersions  in  the  forced  inclinations  from 
the  mean  forced  inclination  obtained  from  the  least  squares  fit  .  .    51 

3.6:  Comparison  of  the  forced  node  obtained  from  the  dynamical  theory 

(bold  line)  and  from  numerical  integration  (points)  of  9  //m 
diameter  Eos-type  dust  particles  as  they  spiral  towards  the  Sun 
from  1.8  AU.  Error  bars  are  the  dispersions  in  the  forced  nodes 
from  the  mean  forced  node  obtained  from  the  least  squares  fit   .  .    52 

3.7:               Comparison  of  the  forced  eccentricity  obtained  from  the  dynamical 
theory  (bold  line)  and  from  numerical  integration  (points)  of  9  /.im 
diameter  Eos-type  dust  particles  as  they  spiral  towards  the  Sun 
from  1.8  AU.  Error  bars  are  the  dispersions  in  the  forced 
eccentricities  from  the  mean  forced  eccentricity  obtained  from  the 
least  squares  fit    53 

3.8:               Comparison  of  the  forced  longitude  of  pericenter  (enforced)  obtained 
from  the  dynamical  theory  (bold  line)  and  from  numerical 
integration  (points)  of  9  /im  diameter  Eos-type  dust  particles  as 
they  spiral  towards  the  Sun  from  1.8  AU.  Error  bars  are  the 
dispersions  in  the  forced  pericenters  from  the  mean  forced 
pericenter  obtained  from  the  least  squares  fit    54 

3.9:  Forced  and  proper  eccentricities  as  functions  of  the  semimajor  axis. 

The  two  straight  reference  lines  have  slopes  of  1.25 55 

3.10:             Variation  in  the  osculating  inclination  vectors  of  9  //m  diameter  Eos 
particles  as  they  spiral  towards  the  Sun.  Each  panel  shows  the 
vectors  at  various  mean  distances  from  the  Sun 56 


vm 


3.11:             Variation  in  the  osculating  eccentricity  vectors  of  9  //m  diameter 
Eos  particles  as  they  spiral  towards  the  Sun,  where  lo  is  the 
longitude  of  pericenter.  Each  panel  shows  the  vectors  at  various 
mean  distances  from  the  Sun 57 

3.12:  Proper  inclination  as  a  function  of  semimajor  axis  of  9  ^l^x\. 

diameter  Eos  dust  particles  as  they  spiral  towards  the  Sun.  Error 
bars  are  the  dispersions  in  the  radii  from  the  mean  radius  obtained 
from  the  least  squares  fit 58 

3.13:             Comparison  of  the  forced  inclination  obtained  from  the  dynamical 
theory  (bold  line)  and  from  numerical  integradon  (points)  of  9  [ivsx 
diameter  Eos  dust  particles  as  they  spiral  towards  the  Sun.  Error 
bars  are  the  dispersions  in  the  forced  inclinations  from  the  mean 
forced  inclination  obtained  from  the  least  squares  fit 59 

3.14:             Comparison  of  the  forced  node  obtained  from  the  dynamical  theory 
(bold  line)  and  from  numerical  integration  (points)  of  9  /<m 
diameter  Eos  dust  particles  as  they  spiral  towards  the  Sun.  Error 
bars  are  the  dispersions  in  the  forced  nodes  from  the  mean  forced 
node  obtained  from  the  least  squares  fit 60 

3.15:             Comparison  of  the  forced  eccentricity  obtained  from  the  dynamical 
theory  (bold  line)  and  from  numerical  integration  (points)  of  9  //m 
diameter  Eos  dust  particles  as  they  spiral  towards  the  Sun.  Error 
bars  are  the  dispersions  in  the  forced  eccentricities  from  the  mean 
forced  eccentricity  obtained  from  the  least  squares  fit 61 

3.16:  Comparison  of  the  forced  longitude  of  pericenter  (unforced)  obtained 

from  the  dynamical  theory  (bold  line)  and  from  numerical  integration 
(points)  of  9  /im  diameter  Eos  dust  particles  as  they  spiral  towards 
the  Sun.  Error  bars  are  the  dispersions  in  the  forced  pericenters  from 
the  mean  forced  pericenter  obtained  from  the  least  squares  fit    .  .    62 


IX 


3.17:  Comparison  of  the  forced  eccentricity  obtained  from  the  dynamical 

theory  (bold  line)  and  from  numerical  integration  (points)  of  9  //m 
diameter  Eos  dust  particles  as  they  spiral  towards  the  Sun.  Error 
bars  are  the  dispersions  in  the  forced  eccentricities  from  the  mean 
forced  eccentricity  obtained  from  the  least  squares  fit.  The  masses 
of  Jupiter  and  Saturn  have  been  reduced  to  0.1  of  their  real  values     63 

3.18:             Variation  in  the  osculating  inclination  vectors  of  9  //m  diameter  Eos 
particles,  in  the  real  solar  system,  as  they  spiral  towards  the  Sun. 
Each  panel  shows  the  vectors  at  various  mean  distances  from  the 
Sun 64 

3.19:             Variation  in  the  osculating  eccentricity  vectors  of  9  i.im  diameter 
Eos  particles,  in  the  real  solar  system,  as  they  spiral  towards  the 
Sun,  where  u;  is  the  longitude  of  pericenter.  Each  panel  shows  the 
vectors  at  various  mean  distances  from  the  Sun 65 

3.20:  Forced  inclination  as  a  function  of  semimajor  axis  for  9  /im 

diameter  Eos  particles  in  1983.  These  are  the  results  of  integrating 
waves  of  particles  from  the  past  such  that  they  end  up  at  various 
distances  from  the  Sun  in  1983 66 

3.21:             Forced  node  as  a  function  of  semimajor  axis  for  9  fim  diameter  Eos 
particles  in  1983.  These  are  the  results  of  integrating  waves  of 
particles  from  the  past  such  that  they  end  up  at  various  distances 
from  the  Sun  in  1983 67 

3.22:  A  vertical  cross-section  of  a  model  dust  band  produced  from  Eos 

particles  alone 68 

3.23:             Forced  inclination  as  a  function  of  semimajor  axis  for  4  different 
systems.  4  G.P.  stands  for  four  giant  planets,  Jupiter,  Saturn, 
Uranus,  and  Neptune.  V,  E,  M,  stand  for  Venus,  Earth,  and  Mars, 
respectively 69 

3.24:             Forced  inclinations  as  a  function  of  semimajor  axis  for  4,  9,  and  14 
/im  diameter  Eos  particles  in  1983.  These  are  the  results  of 
integrating  waves  of  particles  from  the  past  such  that  they  end  up  at 
various  distances  from  the  Sun  in  1983    70 


3.25:             Forced  nodes  as  a  function  of  semimajor  axis  for  4,  9,  and  14  fim 
diameter  Eos  particles  in  1983.  These  are  the  results  of  integrating 
waves  of  particles  from  the  past  such  that  they  end  up  at  various 
distances  from  the  Sun  in  1983 71 

3.26:             Distribution  in  inclination  space  of  particles  with  4  to  9  ^m 
diameters  from  Eos,  Themis,  and  Koronis  families  in  the 
Earth-crossing  region.  The  outer  circle  contains  Eos  particles  while 
the  inner  group  contains  Themis  and  Koronis  particles 72 

3.27:  Seasonal  variation  of  the  relative  numbers  of  asteroidal  particles  that 

will  encounter  the  Earth.  These  are  particles  from  the  Eos,  Themis, 
and  Koronis  families  with  diameters  in  the  range  4  to  9  //m    ...    73 

3.28:  Relative  numbers  of  particles  that  will  encounter  the  Earth  in 

November.  These  are  particles  with  diameters  in  the  range  4  to  9 
fj,m  originating  from  Themis  (bold  line),  Eos  (bold  line),  and  all 
other  asteroids  (thin  line) 74 

3.29:  Variation  in  eccentricity  and  inclination  of  9  /fm  diameter  Eos 

particles  as  they  approach  and  pass  the  Earth.  The  numbers  at  the 
upper  right  hand  corners  of  the  panels  indicate  the  time,  for 
example,  1983-3100  means  3100  years  before  1983 75 

3.30:  Variation  in  the  semimajor  axes  of  fifty  9  /.im  diameter  Eos 

particles  as  they  approach  and  pass  the  Earth.  A  few  of  them  are 
trapped  in  mean  motion  resonances  with  the  Earth.  While  they  are 
trapped,  their  semimajor  axes  remain  constant  and  produce 
horizontal  lines  in  the  diagram 76 

4.1:  Variation  in  the  osculating  eccentricity  vectors  of  9  /<m  diameter 

Encke-type  dust  particles  as  they  spiral  towards  the  Sun,  where  u  is 
the  longitude  of  pericenter.  The  time  interval  between  consecutive 
panels  is  600  years 90 

4.2:  Variation  in  the  osculating  inclination  vectors  of  9  //m  diameter 

Encke-type  dust  particles  as  they  spiral  towards  the  Sun.  The  time 
interval  between  consecutive  panels  is  600  years 91 


XI 


4.3:               The  variation  in  (a,e)  for  every  1000  years  of  one  wave  of  9  //m 
diameter  Encke-type  dust  particles  as  they  spiral  towards  the  Sun. 
A  few  of  them  have  their  semimajor  axes  increased  up  to  4  AU  due 
to  the  process  of  release  from  the  comet 92 

4.4:  Relative  rate  of  change  of  semimajor  axis  (in  arbitrary  units)  as  a 

function  of  eccentricity  due  to  the  action  of  drag  force 93 

4.5:  Distribution  in  (a,e)  space  for  9  /.im  diameter  Encke-type  dust 

particles  in  1983  from  the  second  model.  These  are  the  results  of 
integrating  particles  released  continuously  (once  every  300  years) 
starting  from  the  year  (1983-5,400) 94 

4.6:  Distribution  in  (a,I)  space  for  9  /im  diameter  Encke-type  dust 

particles  in  1983  from  the  second  model.  These  are  the  results  of 
integrating  particles  released  continuously  (once  every  300  years) 
starting  from  the  year  (1983-5,400) 95 

4.7:               Longitude  of  the  ascending  node  (fl)  vs.  longitude  of  pericenter  (a;) 
for  9  //m  diameter  Encke-type  dust  particles  in  1983  from  the 
second  model.  These  are  the  results  of  integrating  particles  released 
continuously  (once  every  300  years)  starting  from  the  year 
(1983-5,400) 96 

4.8:  Distributions  of  the  osculating  eccentricity  vectors  in  different 

semimajor  axis  ranges  for  9  fim  diameter  Encke-type  dust  particles 
in  1983,  from  the  second  model.  The  longitude  of  pericenter  is 
represented  by  u.  These  are  the  results  of  integrating  particles 
released  continuously  (once  every  300  years)  starting  from  the  year 
(1983-5,400) 97 

4.9:  Distributions  of  the  osculating  inclination  vectors  in  different 

semimajor  axis  ranges  for  9  /tm  diameter  Encke-type  dust  particles 
in  1983,  from  the  second  model.  These  are  the  resuhs  of 
integrating  particles  released  continuously  (once  every  300  years) 
starting  from  the  year  (1983-5,400) 98 


xu 


4.10:  Distribution  of  the  osculating  inclination  vectors  in  different 

semimajor  axis  ranges  for  9  jim  diameter  Encke-type  dust  particles 
in  1983,  from  the  first  model.  These  are  the  results  of  integrating 
particles  released  continuously  (once  every  300  years)  starting  from 
the  year  (1983-5,400) 99 

4.11:  Distribution  of  the  inclinations  in  1983  of  9  fim  diameter 

Encke-type  dust  particles  from  the  first  model 100 

4.12:  Distribution  of  the  inclinations  in  1983  of  9  fim  diameter 

Encke-type  dust  particles  from  the  second  model 101 

4.13:  Relative  number  of  orbits  per  unit  semimajor  axis  for  9  fim 

diameter  Encke-type  dust  particles    102 

4.14:  Relative  number  of  orbits  per  unit  semimajor  axis  for  9  yum 

diameter  Eos  dust  particles 103 

4.15:             Forced  eccentricity  as  a  function  of  semimajor  axis  for  9  ^m 
diameter  Encke-type  dust  particles  in  1983.  Error  bars  are  the 
dispersions  in  the  forced  eccentricities  from  the  mean  forced 
eccentricity  obtained  from  the  least  squares  fit    104 

4.16:  Proper  eccentricity  as  a  function  of  semimajor  axis  for  9  /zm  diameter 

Encke-type  dust  particles  in  1983.  Error  bars  are  the  dispersions  in 
the  radii  from  the  mean  radius  obtained  from  the  least  squares  fit  .105 

5.1:  Observed  shape  of  the  zodiacal  cloud  obtained  from  IRAS  data  in 

the  25fim  waveband  with  an  90°  elongation  angle  (SOP406)  .  .  .118 

5.2:               Observed  latitudes  of  the  peak  fluxes  of  the  zodiacal  cloud  around 
the  sky  obtained  from  IRAS  data  in  the  25//m  waveband  with  an 
90°  elongation  angle  (SOP406) 119 

5.3:               Comparison  between  IRAS'  observations  (bold  line)  and  the 
all-asteroid  model  (thin  line)  at  SOP406,  90°  elongation  angle, 
leading  direction 120 


XIU 


5.4:               Comparison  between  IRAS'  observations  (bold  line)  and  the 
all-asteroid  model  (thin  line)  at  SOP392,  97°  elongation  angle, 
leading  direction 121 

5.5:               Latitudes  of  the  peak  fluxes  around  the  sky  from  the  all-asteroid 
model  (thin  curves)  and  the  best  fit  observational  data  (bold 
curves)    122 

5.6:  The  zodiacal  cloud  at  SOP392;  the  bold  lines  are  the  IRAS 

observations.  The  thin  lines  are  from  the  all-asteroid  model  (top), 
the  first  all  comet  model  (central),  and  a  combination  of  asteroid 
and  comet  model  (bottom) 123 

5.7:  The  zodiacal  cloud  at  SOP406;  the  bold  lines  are  the  IRAS 

observations.  The  thin  lines  are  from  the  all-asteroid  model  (top), 
the  first  all  comet  model  (central),  and  a  combination  of  asteroid 
and  comet  model  (bottom) 124 

5.8:  Latitudes  of  the  peak  fluxes  around  the  sky  for  the  27%  asteroid 

plus  73%  comet  model 125 

5.9:  Comparison  between  the  second  all-comet  model  (thin  lines)  and 

the  observations  (bold  lines)  at  two  different  SOP  positions  .  .  .  .126 

5.10:  Inclination  distribution  of  Jupiter  Family  short  period  comets  .  .    127 

5.11:  Distribution  in  (e,  I)  space  of  Jupiter  Family  short  period  comets 

which  have  their  apocenters  inside  the  orbit  of  Jupiter 128 


XIV 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 

of  the  University  of  Florida  in  Partial  Fulfillment  of  the 

Requirements  for  the  Degree  of  Doctor  of  Philosophy. 

DYNAMICAL  EVOLUTION  OF  ASTEROID AL  AND  COMETARY 
PARTICLES  AND  THEIR  CONTRIBUTION  TO  THE  ZODIACAL  CLOUD 

By 

Jer-Chyi  Liou 
December,  1993 

Chairman:   Dr.    Stanley  F.  Dermott 
Major  Department:   Astronomy 

Asteroids  and  comets  are  the  two  major  sources  of  the  interplanetary  dust 
particles  that  populate  the  zodiacal  cloud.  But  how  the  particles  evolve  and  in 
exactly  what  proportion  asteroids  and  comets  contribute  to  this  cloud  are  questions 
that  have  never  been  convincingly  answered.  In  this  dissertation  we  try  to  answer 
these  quesrions  by  studying  the  dynamical  evolution  of  asteroidal  and  cometary 
particles  and  constructing  a  physical  model  to  compare  with  observations  from  the 
Infrared  Astronomical  Satellite  (IRAS). 

We  first  study  the  orbital  evolution  of  asteroidal  and  cometary  particles  under 
the  influence  of  planetary  perturbations,  radiation  pressure,  Poynting-Robertson  drag, 
and  corpuscular  drag.  Analytically,  we  complete  the  development  of  a  "dynamical" 

XV 


secular  perturbation  theory  for  asteroidal  particles.  We  use  a  numerical  integrator, 
RADAU,  on  an  IBM  ES/9000  supercomputer  to  study  the  orbital  evolution  of 
thousands  of  particles,  from  both  asteroids  and  comets,  systematically,  from  their 
origins  to  a  heliocentric  distance  of  around  0.1  AU.  We  then  analyze  the  results  and 
apply  the  distributions  of  the  orbital  elements  from  the  dynamical  study  to  a  three- 
dimensional  model,  SIMUL,  which  generates  the  observed  flux  from  a  given  number 
of  orbits  with  a  known  distribution  of  orbital  elements,  and  compare  the  results  with 
the  IRAS  observations.  We  conclude  that  it  is  possible  to  construct  a  zodiacal  cloud 
model  to  fit  the  shape  of  the  observed  cloud  by  a  combination  of  approximately  27% 
asteroidal  particles  and  73%  cometary  particles. 
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CHAPTER  1 
INTRODUCTION 


The  Zodiacal  Cloud 


The  zodiacal  cloud  is  a  feature  of  the  night  sky  that  can  be  seen  by  the  naked 
eye  just  before  sunrise  and  just  after  sunset.  Its  brightness  increases  towards  the 
ecliptic  plane  and  is  roughly  symmetric  with  respect  to  that  plane.  Since  the  all  sky 
survey  in  1983  by  the  Infrared  Astronomical  Satellite  (IRAS),  the  structure  of  the 
zodiacal  cloud  in  four  infrared  wavebands  (at  12,  25,  60,  and  100  ^m)  has  been 
known  in  great  detail,  including  the  overall  shape  of  the  cloud  and  the  latitudes 
of  the  peak  fluxes  around  the  sky  (e.g.  Hauser  et  al.,  1985).  The  infrared  flux  is 
mainly  due  to  thermal  emission  by  micron-sized  or  larger  dust  particles.  What  are 
the  origins  of  these  particles?  How  do  they  form  and  how  do  they  evolve?  What  are 
the  global  and  fine-scale  structures  of  the  cloud?  Is  there  any  time  variation  in  the 
cloud  structure?  The  observational  data  from  IRAS  has  provided  us  with  a  detailed 
look  at  the  structure  of  the  cloud,  not  only  increasing  our  understanding  of  the  cloud 
but  also  setting  critical  constraints  on  the  models  which  will  be  used  to  explain  the 
formation  and  structure  of  the  cloud. 

Previous  models  of  the  structure  of  the  zodiacal  cloud  are  based  on  some  assumed 
number  density  functions  which  depend  on  heliocentric  distance  and  ecliptic  latitude 
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(e.g.  Giese  et  al.,  1985,  1986,  Kneissel  and  Mann,  1991).  These  models  do  not 
consider  the  formation,  source,  and  evolution  of  the  particles;  they  are  restricted  to 
a  description  of  the  observations  rather  than  an  explanation  of  the  physical  structure 
and  evolution  of  the  cloud.  When  compared  with  IRAS  observations,  none  of  these 
models  are  capable  of  explaining  several  key  features  in  the  data. 

The  dynamical  evolution  of  individual  micron-sized  dust  particles  in  the  solar 
system  has  been  studied  previously  (e.g.  Gustafson  et  al.,  1987,  Jackson  and 
Zook,  1992).  A  more  recent  study  by  Dermott  et  al.  (1992)  has  shown  that  it 
is  possible  to  study  the  orbital  evolution  of  dust  particles  coming  from  different 
sources  systematically  by  incorporating  the  concept  of  the  so-called  proper  and 
forced  elements  from  the  secular  perturbation  theory  (to  be  explained  in  Chapter 
3).  However,  the  work  of  Dermott  et  al.  needs  some  improvement;  in  particular, 
it  did  not,  for  example,  include  the  effects  of  solar  wind  corpuscular  drag  and  the 
changes  in  the  orbital  elements  of  dust  particles  due  to  the  process  of  release.  As 
we  will  show,  these  are  effects  that  significantly  affect  the  orbital  evolution  of  dust 
particles. 

In  summary,  our  aim  is  to  study  the  orbital  evolution  of  dust  particles  coming 
from  different  sources  systematically,  and  including  all  the  important  dynamical 
effects.  We  find  the  structure  of  a  model  zodiacal  cloud  based  on  this  kind  of 
dynamical  study  and  compare  our  results  with  the  IRAS  observations. 
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Objective  of  this  Dissertation 

In  this  dissertation,  we  consider  the  dynamical  evolution  of  dust  particles, 
systematically,  from  source  to  sink,  as  described  by  Dermott  et  al.  (1992)  and 
use  the  results  from  the  dynamical  study  to  construct  a  zodiacal  cloud  model.  We 
know  that  asteroids  and  comets  are  the  two  major  sources  of  the  interplanetary  dust 
particles  that  populate  the  zodiacal  cloud.  But  in  what  way  these  particles  evolve 
and  in  exactly  what  percentage  they  contribute  to  this  cloud  are  questions  that  have 
not  been  answered  before.  These  are  the  questions  that  we  would  like  to  answer. 

In  Chapter  2  we  examine  the  classical  secular  perturbation  theory  which  describes 
the  orbital  evolution  of  planets  due  to  their  mutual  gravitational  perturbations  and 
the  orbital  evolution  of  asteroidal  particles  due  to  the  gravitational  perturbations  of 
the  planets.  We  first  study  the  evolution  of  the  solar  system  numerically  and  find  the 
eigenfrequencies,  amplitudes,  and  phases  associated  with  the  planets  in  the  system 
that  describe  the  variations  of  their  eccentricities  and  inclinations  and  compare  them 
with  predictions  from  secular  perturbation  theory.  We  also  examine  numerically  the 
orbital  evolution  of  asteroidal  particles  in  a  (Sun- Jupiter-Saturn)  system  and,  again, 
compare  the  results  with  those  from  secular  perturbation  theory. 

In  Chapter  3  we  revise  the  classical  secular  perturbation  theory  to  include  the 
effect  of  drag.  The  idea  that  we  use  was  first  developed  by  Dermott  et  al.  in  1992.  We 
complete  the  theory  and  apply  this  "dynamical"  theory  to  the  evolution  of  asteroidal 
dust  particles  starting  from  several  different  heliocentric  distances  until  they  reach 
the  inner  part  of  the  solar  system,  we  then  compare  these  analytical  results  with  the 
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numerical  results.  We  next  integrate  numerically  the  orbits  of  thousands  of  4,  9, 
and  14  /.im  diameter  asteroidal  dust  particles  systematically  from  their  origin  until 
they  reach  about  0.1  AU  from  the  Sun  in  the  real  solar  system.  All  the  planetary 
perturbations  (except  those  of  Mercury  and  Pluto),  radiation  pressure,  Poynting- 
Robertson  light  drag,  and  solar  wind  corpuscular  drag  effects  are  included  in  our 
calculations.  We  analyze  the  data  and  obtain  the  distributions  of  their  orbital  elements 
as  functions  of  semimajor  axis,  and  also  examine  the  orbits  of  these  dust  particles 
when  they  approach  the  Earth  to  gain  information  useful  for  future  Earth-orbiting 
dust  collection  facilities. 

In  Chapter  4  we  study  the  evolution  of  Encke-type  dust  particles  using  procedures 
similar  to  those  in  Chapter  3.  This  is  the  first  time  that  this  has  been  attempted  for 
particles  in  cometary  type  orbits.  We  construct  two  models  in  which  the  dust  particles 
are  given  different  initial  inclination  distributions.  We  follow  their  orbital  evolution 
numerically  as  they  spiral  towards  the  Sun.  The  data  are  analyzed  in  a  similar  way 
to  those  in  Chapter  3  in  order  to  obtain  the  distributions  of  their  orbital  elements 
as  functions  of  semimajor  axis.  The  difference  in  the  initial  inclinations  turns  out 
to  be  a  key  element  in  constructing  the  zodiacal  cloud  model.  We  also  show  that, 
because  of  their  high  eccentricities,  comets  cannot  be  the  source  of  any  of  the  solar 
system  dust  bands  discovered  by  IRAS. 

We  use  the  results  from  Chapters  3  and  4  to  construct  several  physical  zodiacal 
cloud  models  (Chapter  5).  A  three-dimensional  model,  SIMUL,  is  used  to  calculate 
the  flux  from  the  zodiacal  cloud  model.    We  compare  the  shape  of  the  zodiacal 
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cloud  and  latitudes  of  the  peak  fluxes  from  different  zodiacal  cloud  models  with  the 
observations  from  IRAS.  Finally  we  show  that  by  combining  cometary  dust  particles 
and  asteroidal  dust  particles,  we  can  construct  a  model  that  describes  the  observed 
shape  of  the  zodiacal  cloud.  We  also  discuss  alternative  models  and  some  unsolved 
problems. 


CHAPTER  2 
SECULAR  PERTURBATION  THEORY 

In  this  chapter  we  will  first  summarize  the  theory  of  secular  perturbation,  which 
is  the  classical  method  to  study  the  long-term  orbital  evolution  of  planets  as  well  as 
minor  planets  in  the  solar  system.  We  apply  this  theory  to  a  (Sun-Jupiter-Saturn) 
system  to  find  the  eigenfrequencies,  amplitudes,  and  phases  associated  with  the 
planets  in  the  system  and  compare  the  results  with  those  from  the  direct  numerical 
integration.  We  then  use  the  results  from  the  numerical  integration  to  calculate  the 
initial  forced  elements  for  asteroidal  type  particles  at  2.6  AU  and  follow  their  orbital 
evolution  in  this  (Sun-Jupiter-Satum)  system  both  numerically  and  analytically  and 
compare  the  results. 

Low-order  Secular  Perturbation  Theory 

The  Solar  System 

The  classical  way  to  study  the  long  term  rates  of  change  of  orbital  elements, 
eccentricities  (e),  longitudes  of  pericenter  (cc),  inclinations  (I),  and  longitudes  of 
ascending  node  (17)  of  small  bodies  in  the  solar  system  without  any  mean  motion 
resonance  involved  is  by  using  secular  perturbation  theory  (Brouwer  and  Clemence, 
1961,  Dermott  et  al,  1986).  It  is  based  on  the  assumption  that  the  eccentricities  and 
inclinations  of  all  objects  involved  are  small.  The  long  term  variations  in  e,  to,  I, 
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and  n  of  a  small  body  are  largely  determined  by  the  secular  terms  in  its  disturbing 
function  (to  be  defined  later).  These  are  terms  independent  of  the  mean  longitudes 
of  any  objects  involved. 

Considering  the  solar  system  with  the  Sun  at  the  center,  the  Lagrange  equations 
which  describe  the  rates  of  change  of  the  orbital  elements  of  any  given  planet  are 

da  __    2  dR 
dt       na  d\ 


de       -VT^^,  / -,dR      VT^^dR 

17  = ^ (^  -  Vl  -  e-j-TTT- ^ — ^— 

dt  na-e  aX  nwe     ow 

dl  _      -tanU     .OR       dR 1  OR 

H  ~  na2v/r^   5A  ^  5^^  ~  na2yr3^sin7M 

dn  _  1  dR  (2.1) 

dt        naVl  -e^sin/  dl 

dw  _  VT^^dR  tani/      dR 

dt  na-e     de       na'^Vl  -  e-  dl 

de_-2dR      \/r^?(l  -  VT^^)  dR  tan  ^7      dR 

dt        na  da  na"e  de       na'^\/l  —  e'^  dl 


(2.2) 


where 

n"a^  =  G{M  +  m) 

A  =  e  +       ndt 

and  G,  M,  m,  n,  R  are  the  universal  gravitational  constant,  mass  of  the  sun,  mass  of 
the  given  planet,  the  mean  motion  of  the  given  planet,  and  the  disturbing  function, 
respectively.  The  disturbing  function  is  the  potential  acting  on  this  given  planet  due 
to  the  existence  of  other  planets  in  the  system.  The  semimajor  axis,  eccentricity, 
inclination,  longitude  of  the  ascending  node,  longitude  of  pericenter,  and  mean 
longitude,  are  represented  by  a,  e,  I,  n,  ta,  and  A  ,  respectively. 
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In  order  to  avoid  singularities  that  arise  in  Equation  2. 1  for  small  eccentricities 
and  inclinations,  we  define  four  quantities  {h,k,p,q)  to  replace  e,  cc?,  I,  and  Q: 

<  k  =  e  cos  w 

p  =  tan  /  sin  $7  ,^  ^, 

fro  2.4 

q  =  tan  I  cosSZ  . 

If  the  eccentricities  and  inclinations  of  the  planets  are  small,  we  can  ignore  the 
second  and  higher  order  terms  in  e  and  I  in  Equation  2. 1  and  the  resultant  Lagrange 
equations  involving  e,  w,  I,  and  fl  become 

n,a~,  dk, 

TZjOj  dhj 

and 

1     dRj 

P^  = 27r^ 

1    dRj  (2-^^ 

iija--  upj 
By  using  similar  assumptions,  the  disturbing  function  to  the  second  order  in  e 
and  I  experienced  by  the  j'*'  planet  has  the  following  form  (Dermott  and  Nicholson, 
1986): 

2    2  V^  "^i  f  1  -     ;(l)/'2         72\         ^  -     .(2)  ,  A 

2    2  V^  '^J   f  1  -       (1)  1 

(2.7) 


or 


(2.8) 


and 


where  j  is  not  equal  to  /  and 

J  aijuj        if  Gj  >  Qj     {internal  perturber) 
■''       ]^  cij/ai       if  dj  <  tti     {external  perturber) 

(  1  if  a  J  >  at     {internal    perturber) 

■'^        [  aj/ui        if  aj  <  Ut     {external    perturber) 

1         ^ 

^      ^=l  ^^  (2.11) 

1  nii  _      10) 

^       t=i  ^'^  (2.12) 

R      -   ^        "^^         -     a(1) 

and  63  A,  63"^!,  are  the  so-called  Laplace  coefficients  (e.g.  Brouwer  and  Clemence, 
1961).  Mercury  and  Pluto  are  not  included  in  our  calculation  due  to  their  small 
masses. 

The  solutions  of  hj,  kj,  pj,  qj  are 

7 
hj  =  Y^  ejism{gjt  +  (3t) 

^  'f  (2.13) 

^]  =  Z  ^ji  cos{git  +  I3i) 

and 

6 

P]  =  12  Ijism{ftt  +  7j) 

'f  (2.14) 

9j  =   E  IjtCOs{fjt  +  -ft) 

1=1 
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where  gi  and  fi  are  the  eigenfrequencies  of  the  matrices  A  and  B  whose  elements 
are  defined  by  Ajj,  Aji,  Bjj,  and  Bji  above.  The  amplitudes,  Cji,  Iji,  and  phases,  /?i,  71 
are  determined  from  the  initial  conditions  of  the  planets. 

Equations  2.13  and  2.14  are  the  results  from  the  classical  low-order  secular  per- 
turbation theory.  These  two  equations  describe  analytically  the  respective  variations 
in  the  eccentricity,  inclination,  longitude  of  pericenter,  and  longitude  of  ascending 
node  of  the  j''^  planet.  They  are  based  on  the  assumption  that  the  eccentricities  and 
inclinations  of  all  the  objects  involved  are  small.  When  there  is  no  mean  motion  res- 
onance phenomenon  involved,  this  theory  describes  fairly  accurately  the  long-term 
variations  in  e,  tu,  I,  and  fl  of  the  planets. 
Small  Bodies  in  the  Solar  System 

For  small  bodies  moving  under  the  gravitational  perturbations  of  planets,  we 
can  follow  similar  procedures  to  find  the  long-term  rates  of  change  in  their  {h,k), 
and  (p,q).    If  we  use  A  to  represent  the  nodal  precession  rate  of  the  small  body 

(subscript  s),  i.e., 

7 

n  s—^  ITT-i  .(1) 

^  =  lEl7^^^^^^V2  (2.15) 

where  n,  mi,  and  M  are  the  mean  motion  of  the  small  body,  the  mass  of  the  i'^  planet, 
and  the  mass  of  the  Sun,  respectively,  then  the  solutions  for  {h,k)  and  {p,q)  are 

/i  =  epsin(/l< +  /?) +  /io 

k  =  tp  cos[  At +  0)  +  kQ  ^      ^ 

p  =  Ip  sm{- At  +  7)  +  po  ^2  17) 

q  =  IpCOs{-At  +  j)  +  qo 


{ 


where 
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7 

sm{gtt  +  /3t) 

A  —  Qi 

7 
^^0  =  -  y^  -r—^ —  cos{g,t  +  0i) 

6 


(2.18) 


1=1 

6 


(2.19) 


50  =  J]  /I   I '  f  cos(/8i  +  7j) 
,=1  ^  +  •^' 

n  J^mi       _      (2) 

/=i 
n^^mi       _      (1) 

1=1 


(2.20) 


and  g,-,/i,  eii,  In,  and  /9,-,  7/  are  the  eigenfrequencies,  amplitudes,  and  phases  associated 
with  the  planets  of  the  solar  system  from  the  previous  section. 

In  the  previous  equations,  ep,  Ip,  fS  and  7  can  be  determined  from  the  initial 
conditions  of  the  small  body.  The  forced  eccentricity  and  inclination  are  defined  as 


^f  =  yJPo  +  % 
The  facts  that  {h,k),  and  ip,q)  are  de-coupled  from  each  other  and,  for  each  of 

them,  we  can  use  the  proper  and  forced  elements  (their  characteristics  under  the 

influence  of  drag  force  will  be  described  in  Chapter  3)  instead  of  the  osculating 

elements  to  describe  their  evolution  make  it  possible  to  study  dust  particles  coming 

from  the  same  source  systematically  (Dermott  et  al.,  1985). 

The  precession  term.  A,  depends  on  the  masses  of  the  planets  if  the  semimajor 

axis  of  the  minor  planet  is  given.  This  means  the  orbit  of  the  minor  planet  precesses 
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due  to  the  existence  of  planets.  If  all  the  planets  have  circular  orbits  and  the  same 
inclinations,  then  there  will  be  no  forced  eccentricity  and  longitude  of  pericenter  (due 
to  zero  eccentricities  of  planets)  and  the  orbit  of  the  minor  planet  will  precess  along 
the  common  plane  of  the  planets  (due  to  the  same  inclinations  of  planets).  However, 
this  is  not  the  case  in  the  real  solar  system.  The  forced  components  in  Equations 
2.16  and  2.17  will  be  present  due  to  the  eccentricities  and  different  inclinations  of 
the  planets'  orbits.  The  forced  components  introduce  a  plane  of  symmetry  (different 
from  any  of  the  planets'  orbital  planes),  with  respect  to  which  the  orbit  of  the  minor 
planet  precesses,  along  with  the  displacement  of  the  pattern  of  the  orbit  of  the  minor 
planet  away  from  the  Sun  (Dermott  et  al.,  1985). 

The  equations  of  (h,k)  and  (p,q)  describe  the  time  variation  in  orbital  elements  of 
small  bodies  in  the  solar  system.  This  is  a  result  from  the  low-order  theory.  When 
the  eccentricity  or  inclination  of  the  small  body  becomes  larger,  the  low-order  theory 
can  no  longer  be  applied  to  study  the  evolution  of  this  minor  planet's  orbit.  In  this 
case,  direct  numerical  integration  is  the  only  way  to  study  its  orbital  evolution.  In  the 
next  section  we  will  test  this  low-order  theory  on  the  real  solar  system  and  compare 
the  results  with  those  from  direct  numerical  integration. 

Numerical  Approach 

After  the  invention  of  high  speed  computers,  astronomers  found  another  powerful 
tool  to  study  celestial  mechanics:  numerical  integration.  Many  effects  that  can 
not  be  handled  easily  by  analytical  theory  can  be  included  in  the  direct  numerical 
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integration.  For  example,  if  we  include  the  gravitational  forces  due  to  the  Sun  and 

planets  and  the  radiation  pressure  force  as  well  as  the  drag  terms  in  the  equation 
of  motion  of  a  dust  particle  and  numerically  integrate  the  equation,  passage  through 
mean  motion  resonance,  planetary  scattering,  and  trapping  in  mean  motion  resonance 
associated  with  planets  can  be  seen  from  the  numerical  output  of  the  orbit  of  this 
particle.  However,  the  numerical  error  built  up  in  the  process  needs  to  be  considered 
carefully.  Besides,  not  all  problems  can  be  solved  numerically.  In  some  cases,  for 
example,  when  the  number  of  particles  involved  becomes  too  large  or  the  integration 
time  becomes  too  long,  even  the  fastest  machine  in  the  world  is  not  sufficient  to 
obtain  the  desired  results. 

The  numerical  integrator  we  chose  is  RADAU.  It  was  developed  by  Everhart 
(1985).  The  major  computer  we  use  to  run  the  numerical  integration  is  an  IBM 
supercomputer  ES/9000  in  the  Northeast  Regional  Data  Center  (NERDC)  at  the 
University  of  Florida.  IBM  ES/9000  is  equipped  with  a  vector  facility.  For  each 
processor,  it  can  take  up  to  256  objects  and  perform  the  same  operation  at  the  same 
time.  In  order  to  improve  the  performance,  we  have  vectorized  RADAU.  In  addition 
to  IBM  ES/9000,  we  have  also  used  IBM  RS/6000  350,  IBM  RS/6000  340,  IBM 
RS/6000  220,  and  several  Sun  workstations. 

The  numerical  integrations  we  describe  in  this  chapter  include  two  parts.  The 
first  part  is  the  integration  of  the  solar  system  with  only  two  planets,  Jupiter  and 
Saturn,  with  their  real  masses  and  orbital  elements.  This  is  a  simple  three-body 
system.  We  have  also  integrated  the  system  with  the  reduced  masses,  eccentricities. 
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and  inclinations  of  Jupiter  and  Saturn.  Both  numerical  results  were  analyzed  and 
compared  with  analytical  predictions  from  the  secular  perturbation  theory.  The 
second  part  is  the  integration  of  the  orbital  evolution  of  100  asteroidal  type  particles 
in  the  (Sun-Jupiter-Saturn)  system.  Again,  the  comparison  between  analytical  and 
numerical  results  was  made. 

Comparison  between  the  Analytical  and  Numerical  Calculations 

Sun-Jupiter-Satum  System 

We  first  test  the  secular  perturbation  theory  on  the  (Sun- Jupiter-Saturn)  system. 
This  three-body  system  was  integrated  from  1983  backward  for  five  hundred  thousand 
years.  The  variations  in  h,  k,  p,  and  q  of  Jupiter  and  Saturn  are  shown  in  Figure  2.1 
and  Figure  2.2.  We  analyzed  the  data  by  using  the  Fast  Fourier  transform  method 
to  find  the  eigenfrequencies,  amplitudes,  and  phases  of  the  system  from  the  power 
spectra  and  compared  the  results  with  the  analytical  calculation.  The  comparison  is 
shown  in  Table  2.1. 

As  we  can  see  from  the  table,  the  two  results  agree  reasonably  well.  However, 
there  are  certain  discrepancies  for  the  eigenfrenquencies,  gi  and  g2.  The  differences 
are  15.2%  and  16.8%,  respectively.  A  possible  reason  for  this  might  be  the  mean 
motion  resonance  between  Jupiter  and  Saturn.  These  two  planets  are  very  close  to 
5:2  mean  motion  commensurability  (a  third  order  resonance).  The  analytical  third- 
order  resonance  theory  has  not  been  worked  out,  and  its  development  is  beyond  the 
scope  of  this  dissertation. 
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We  can  also  check  this  mean  motion  resonance  effect  by  reducing  the  masses, 
eccentricities,  and  inclinations  of  Jupiter  and  Saturn  and  compare  the  results  again. 
The  comparison  is  in  Table  2.2.  The  masses  of  Jupiter  and  Saturn  are  10%  of 
their  original  masses.  Their  eccentricities  and  inclinations  are  ej  =  0.0004948,  Ij  = 
0.003273  degree,  es  =  0.0005085,  and  Is  =  0.008872  degree.  The  largest  difference 
is  only  3.3  %.  All  the  other  differences  are  less  than  1%.  The  secular  permrbation 
theory  predicts  very  accurately  the  results  from  the  numerical  integrations  in  this  case. 

From  these  comparisons  we  can  conclude  that  the  low-order  secular  permrbation 
theory  does  describe  the  behavior  of  planets  accurately  except  when  a  mean  motion 
resonance  between  planets  is  involved.  We  will  further  examine  how  good  this  theory 
is  when  we  apply  it  to  the  asteroidal-type  particles  in  the  next  section. 

Asteroidal  Particles  in  the  Solar  System 

In  the  calculations  described  in  this  section  we  put  100  asteroidal-type  particles  in 
the  Sun- Jupiter-Saturn  system  and  integrate  their  evolution  for  280,000  years.  Again, 
we  compare  the  results  with  the  analytical  calculation  from  secular  perturbation 
theory. 

The  semimajor  axes  for  the  particles  are  chosen  to  be  2.6  Astronomical  Units 
(AU).  This  is  to  avoid  several  mean  motion  resonance  locations  with  Jupiter  at 
around  3  AU.  We  use  the  "particles  in  a  circle"  method  (Dermott  et  al.,  1985,  1992) 
to  set  up  the  initial  orbital  elements  of  the  particles.  Their  proper  eccentricities  and 
inclinations  are  0.049  and  2.12  degrees,  respectively,  which  are  the  proper  elements 
of  the  Koronis  family.    The  proper  Q,  w,  and  mean  longitudes  of  particles  are 
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randomly  chosen  between  0°  and  360°.  The  initial  forced  elements  at  2.6  AU  are 
calculated  by  using  the  eigenfrequencies,  amplitudes,  and  phases  from  the  numerical 
integration  of  the  (Sun- Jupiter-Saturn)  system  from  the  previous  section. 

The  positions  of  particles  in  the  (p,q)  and  {h,k)  phase  spaces  at  t  =  0  and  t  = 
280,000  years  are  shown  in  Figures  2.3  and  2.4.  Due  to  the  difficulty  in  generating  w 
in  the  diagram,  we  use  u  to  represent  the  longitude  of  pericenter  in  all  diagrams.  The 
fact  that  after  280,000  years  all  the  particles  still  remain  in  pretty  well-defined  circles 
in  both  diagrams  shows  that  the  method  of  using  the  forced  and  proper  elements  to 
describe  the  evolution  of  asteroidal-type  particles  is  a  good  one. 

The  variations  in  the  forced  elements  are  shown  in  Figures  2.5  and  2.6.  The 
dots  are  results  from  the  direct  numerical  integrations  while  the  solid  lines  are  results 
from  the  analytical  calculation  from  the  secular  perturbation  theory.  The  differences 
are  about  0.2°  in  the  forced  inclination  and  less  than  0.01  in  the  forced  eccentricity. 
These  results  are  actually  not  too  bad.  This  shows  that  when  the  particles  are  not  in 
mean  motion  resonance  with  any  planet,  the  secular  perturbation  theory  does  describe 
the  evolution  of  asteroidal-type  particles  to  a  certain  accuracy.  However,  when  small 
particles  are  spiralling  towards  the  Sun  due  to  Poynting-Robertson  and  corpuscular 
drag,  this  classical  secular  perturbation  theory  does  not  apply  any  more.  In  the  next 
chapter  we  use  the  theory  of  Dermott  et  al.  (1992)  to  include  the  effect  of  drag. 
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CHAPTER  3 
ORBITAL  EVOLUTION  OF  ASTEROIDAL  PARTICLES 


Radiation  Pressure,  Poynting-Robertson  Drag,  and  the  Corpuscular  Drag 

The  orbital  evolution  of  micron-sized  dust  particles  in  the  solar  system  is  affected 
by  several  forces  in  addition  to  the  gravitational  forces  due  to  the  Sun  and  planets. 
These  include  radiation  pressure,  Poynting-Robertson  (PR)  drag,  and  solar  wind 
corpuscular  drag. 

The  radiation  pressure  force  is  the  force  exerted  on  a  dust  particle  due  to  the 
interception  of  sunlight.  It  varies  as  the  inverse  square  of  the  distance  from  the  Sun 
and  depends  on  the  properties  of  the  particle  (size,  chemical  composition,  etc.). 

When  a  micron-sized  dust  particle  is  moving  around  the  Sun,  it  absorbs  the 
sunlight  and  reradiates  the  energy  away  isotropically  in  its  own  reference  frame. 
However,  when  viewed  from  the  rest  frame  this  particle  radiates  more  energy  along 
the  same  direction  as  its  velocity  vector.  This  will  cause  a  net  loss  in  its  orbital 
energy  and  angular  momentum.  It  acts  like  a  retarding  force  on  the  particle.  This 
deceleration  drag  causes  the  particle  to  spiral  towards  the  Sun.  It  is  called  the 
Poynting-Robertson  drag  (Poynting,  1903,  Robertson,  1937).  This  effect  is  present 
even  if  the  particle  scatters  the  sunlight  instead  of  absorbing  and  reradiating  the 
radiation. 
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The  solar  wind  corpuscular  drag  is  due  to  the  interaction  between  the  solar  wind 
protons  and  the  dust  particles.  The  physical  interpretation  is  similar  to  that  of  the  PR 
drag.  The  classical  paper  dealing  with  radiation  pressure,  PR  drag,  and  corpuscular 
drag  is  the  paper  by  Burns  et  al.    (1979). 

The  equation  of  motion  of  a  dust  particle  with  geometrical  cross  section  A  and 
mass  m  moving  under  the  influence  of  the  gravitational  forces  of  the  Sun,  planets, 
and  the  radiation  force  from  the  solar  radiation  can  be  written  as 


mv 


GMqUI 


9 


EGMnm 


n=l         '■"  (3.1) 

l-^>-^ 

c  /  c 


+  {SA/c)Qpr 

where  v  and  r  are  the  velocity  and  posidon  vectors  of  the  particle  in  the  heliocentric 
coordinate  system,  while  M„,  Pn  are  the  mass  of  the  n''^  planet  and  the  position 
vector  of  the  particle  with  respect  to  that  planet,  respectively.  S,  c,  and  Qpr  are  the 
solar  energy  flux  density,  the  speed  of  light,  and  the  radiation  pressure  coefficient, 
respectively.  The  particle's  radial  velocity  and  the  unit  vector  of  the  incident  radiation 
are  represented  by  r  and  S,  respectively.  The  non-velocity  dependent  part  of  the  third 
term  at  the  right  hand  side  of  Equation  3.1  is  the  radiation  pressure  force  term,  while 
the  velocity  dependent  parts  are  the  drag  terms.  Qpr  depends  on  the  properties 
(density,  shape,  size,  etc.)  of  the  particle.  It  can  be  calculated  from  the  Mie  theory 
(Bums  et  al.,  1979,  Gustafson,  1994).  Traditionally,  a  dimensionless  quantity,  /3,  is 
introduced  to  specify  the  effect  of  radiation  pressure  and  PR  drag.  It  is  defined  as 

gravitation   jorce 
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The  major  topic  in  this  chapter  is  to  study  the  dynamical  evolution  of  asteroidal 
dust  particles  which  contribute  most  to  the  25  /im  waveband  IRAS  observations.  If 
we  adopt  the  solar  wind  drag  as  being  35%  of  the  PR  drag,  the  major  contributors  to 
the  25  fim  waveband  are  particles  with  f3  equal  to  0.05037,  corresponding  to  9  fim 
in  diameter.  These  results  are  obtained  by  assuming  2.7  gm/cm^  perfectly  spherical 
astronomical  silicate  and  using  Mie  theory  to  integrate  over  the  solar  spectrum 
(Gustafson,  1994).  The  corpuscular  drag  is  35%  that  of  the  PR  drag  according 
to  the  latest  average  solar  wind  data  from  Helios  1  and  2  (Schwenn,  1990). 

When  a  dust  particle  is  released  from  a  large  parent  body,  it  immediately  "sees" 
a  less  massive  Sun,  due  to  the  radiation  pressure.  Because  the  dust  particle  has 
identical  position  and  velocity  vectors  as  its  parent  body  (assuming  it  leaves  with 
zero  relative  velocity),  this  reduction  in  the  centripetal  force  will  change  some  of  its 
orbital  elements  right  away.  Since  radiation  force  is  a  radial  force,  it  has  no  effect 
on  /  and  H  of  the  dust  particle;  only  the  semimajor  axis,  eccentricity,  and  longitude 
of  pericenter  need  to  be  recalculated. 

The  new  semimajor  axis,  an,  and  eccentricity,  en,  are  given  by 

1-/3 


On  —  a 


671   — 


1  - 


1  -2a(3/r 
(l-2a^/r)(l-e2) 


(3.3) 


1 

2 

1 


=  X2  (3.4) 


{i-i3r 

where  a  and  e  are  the  semimajor  axis  and  eccentricity  of  the  parent  body  and  r  is 
the  radial  distance  from  the  sun  at  which  the  dust  particle  is  released.  From  the  first 
equation  we  can  see  that  the  new  semimajor  axis  of  the  dust  particle  is  always  larger 
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than  that  of  its  parent  body.  However,  this  is  not  true  for  the  new  eccentricity.  The 
new  semimajor  axis  and  eccentricity  of  the  dust  particle  also  depend  on  where  it  is 
released  (i.e.,  r).  It  is  not  too  difficult  to  show  that  both  a„  and  e„  have  their  maxima 
at  perihelion  release  and  have  their  minima  at  aphelion  release.  When  the  quantity 
X  in  Equation  3.3  becomes  negative,  the  new  eccentricity  is  the  square  root  of  the 
absolute  value  of  X,  and  there  is  an  180°  shift  in  the  longitude  of  pericenter  of  the 
dust  particle  compared  with  that  of  its  parent  body. 

Static  Theory  and  Dynamical  Theory 

The  classical  secular  perturbation  theory  mentioned  in  Chapter  2  deals  with 
particles  having  constant  semimajor  axes.  This  is  the  so-called  static  theory.  Micron- 
sized  dust  particles  in  the  solar  system  will  spiral  in  towards  the  Sun  due  to  the  drag 
force  in  about  50,000  years  (Wyatt  and  Whipple,  1950,  Burns  et  al,  1979).  Thus,  it 
is  necessary  to  develop  a  new  secular  perturbation  theory  which  includes  the  effect 
of  drag.  We  use  the  theory  developed  in  Dermott  et  al.  1992.  We  test  the  theory 
on  asteroidal  dust  particles  in  this  section. 

In  the  low-order  static  theory,  the  rates  of  change  of  p,  q,  h,  and  k  are 

h  =  Ak  —  ^  Uj  cos{gjt  +  (3j) 


j=l 


N  (3.5) 


k  =  —Ah  +  ^  Vj  s\n{gjt  +  ^j) 

N 
p  =  -Aq  +  Y.  /'j  cos(/ji  +  7j) 


M  (3-6) 


q  =  Ap-Y,  ^ij  s\n{fjt  +  jj) 
3=1 
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where  A,  uj,  j.ij,  gj,  fj,  /3j,  and  7y,  are  defined  in  Chapter  2.  When  there  is  no  drag, 
the  precession  rate,  A,  is  a  constant.  When  there  is  a  drag  force,  a  dust  particle 
spirals  in  towards  the  Sun.  Its  semimajor  axis  becomes  a  function  of  time.  So,  A, 
Uj,  and  i-Lj  (through  their  dependences  on  semimajor  axis)  all  become  functions  of 
time.  The  eigenfrequencies  and  phases,  gj,  fj,  /3j,  and  fj,  are  constants  regardless 
of  the  presence  of  drag  force. 

In  order  to  solve  the  problem,  we  define  a  complex  equation  to  combine  the 
equations  involving  p  and  q  (or  h  and  k),  namely 

z  =  q  +  ip  .  (3.7) 

By  using  z.  Equation  3.6  becomes 

N 
i  =  -iAz  +  i  Y^  fije'^^^^+''^^  .  (3.8) 

In  the  static  theory,  the  solution  of  the  above  equation  is 

N 
z  =  ^e'(-^'+^)  +  y  -J^e^if^'+-'^)  (3.9) 

where  //  and  7  are  the  proper  inclination  and  node,  respectively,  of  the  particle  at 

time  equals  to  zero. 

When  A  and  fij  become  functions  of  time.  Equation  3.8  becomes  a  standard  first 

order  inhomogeneous  differential  equation.  Its  solution  is 

t=tf 


z{t)  =  e- 


C+  e'rdt 

t=0 


(3.10) 
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where 


t=tf 
s=    f  lAdt  (3.11) 

and 

N 
r  =  z^/ije'(^^'+^^)  .  (3.12) 

The  constant  C  in  Equation  3.10  depends  on  the  initial  conditions  of  the  dust  particle. 
Before  we  find  the  solution  of  Equation  3.10,  we  need  to  know  the  rates  of 
change  of  semimajor  axis  and  eccentricity  due  to  PR  drag.   They  are  (Wyatt  and 
Whipple,  1950,  Burns  et  al.,  1979) 

dt)  (l-e2)3/2 

(3.14) 


de\  —5r]QpTe 


dt)       2a2(l-e2)'/2 
where  rj  =  2.53  x  10^^ /ps  for  spherical  particles  in  heliocentric  orbits  with  density 

p,  radius  s  in  (c.g.s.)  units,  and  Qpr  is  the  radiation  pressure  coefficient. 

Now  let  us  take  a  look  at  how  the  drag  force  affects  the  forced  and  proper 

eccentricities.  We  first  consider  a  group  of  dust  particles  from  a  given  family  having 

a  distribution  in  the  (h,k)  phase  space  as  shown  in  Figure  3.1.  The  equation  of  the 

circle  is  given  by 

(k-kcf  +  ih-hS'=elo  (3.15) 

where  he  and  kc  define  the  forced  elements  (i.e.,  the  center  of  the  circle)  at  time  equal 
to  0  and  epo  is  the  initial  proper  eccentricity.  The  total  eccentricity  of  a  given  particle 
in  the  circle  at  time  equal  to  0  is  defined  by  eto-  After  a  small  time  interval,  St,  the 
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center  of  the  circle  moves  to  ni  while  the  radius  of  the  circle  changes  to  epi.  The 
forced  eccentricity  changes  from  efo  to  efi.  The  change  in  the  proper  eccentricity 
(i.e.,  the  radius  of  the  circle)  is  due  to  drag,  while  the  change  in  the  forced  eccentricity 
and  longitude  of  pericenter  is  due  to  both  gravitational  perturbation  and  drag.  From 
Equadon  3.14  we  see  that  the  drag  force  changes  the  total  eccentricity  of  any  given 
particle  by  the  amount 

6e  =  -CietQ  (3.16) 

where  C;  is  a  posidve  constant.  The  new  total  eccentricity  of  the  previous  particle 
in  the  circle,  Cii  then  becomes 

^n  (=  \f^k+hl^  =  (1  -  Ci)e,o   .  (3.17) 

The  equation  of  the  new  circle  after  ^t  becomes 

[(1  -  Ci)k  -  (1  -  Ci)fcJ'  +  [(1  -  Ci)h  -  (1  -  Ci)/ic]2 

(3.18) 
=  [{l-Ci)ej,af  . 

This  means  that  the  drag  force  moves  the  center  of  the  circle  (which  is  defined  by 
the  forced  eccentricity  and  longitude  of  pericenter)  towards  the  origin  by  the  factor 
of  (l-Cy)  and,  at  the  same  time,  makes  the  circle  (defined  by  its  radius,  the  proper 
eccentricity)  shrink  by  the  same  factor.  This  factor,  can  of  course  be  calculated 
analytically  from  Equation  3.14. 

The  overall  picture  of  the  dust  evolution  thus  becomes  very  clear.  The  gravita- 
tional permrbation  changes  the  forced  eccentricities  of  the  dust  particles  while  having 
no  effect  on  their  proper  eccentricities.  The  drag  force  makes  the  dust  particles  spiral 
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towards  the  Sun  and  changes  their  proper  and  forced  eccentricities  the  same  way  it 
changes  their  total  eccentricities. 

By  combining  Equations  3.13  and  3.14  we  can  obtain  the  following  equation: 

de        5e(l  —  e") 

—  =  — ^^ (3  19") 

da       2a{2  +  3e2)  "  ^  "     ^ 

For  asteroidal  dust  particles  we  can  drop  the  e^  terms.  The  slope  in  a  log(a)  vs. 
log(ep)  plot  during  their  evolution  is  then  1.25.  This  will  be  verified  numerically 
in  the  next  section. 

Numerical  Comparisons 

In  this  section  we  first  test  the  dynamical  theory  by  placing  dust  particles  at  1.8 
AU  initially  in  a  (Sun- Jupiter-Saturn)  system.  This  is  to  avoid  the  passage  through 
mean  motion  resonances  associated  with  Jupiter.  Then  we  place  dust  particles 
initially  at  3  AU  and  observe  how  they  are  affected  when  they  do  go  through  the 
Kirkwood  Gaps  (to  be  defined  later). 

Drag  Effect 

We  place  254  Eos-type  dust  particles  initially  at  1.8  AU  to  start  the  numerical 
integration.  Again,  we  use  the  "particles  in  a  circle"  method  (Dermott  et  al.,  1992)  to 
set  up  their  initial  orbital  elements.  The  particles  have  the  same  proper  eccentricities 
and  inclinations  as  those  of  Eos.  Their  initial  forced  elements  are  calculated  using 
the  eigenfrequencies,  amplitudes,  and  phases  from  the  numerical  integration  of  the 
(Sun-Jupiter-Saturn)  system.  We  did  not  include  the  change  in  the  orbital  elements 
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(a,  e,  and  zu)  due  to  the  releasing  process.   The  reason  for  this  is  to  simplify  the 

problem  and  study  the  effect  of  drag  on  the  orbital  evolution  alone. 

The  dust  particles  are  integrated  for  20,000  years.  They  end  up  around  0.5  AU. 
We  plot  their  (p,q)  and  (h,k)  once  every  2,200  years  as  shown  in  Figure  3.2  and 
Figure  3.3.  The  fact  that  they  remain  in  nice  circles  in  both  (p,q)  and  ih,k)  phase 
spaces  during  their  evolution  proves  that  the  method  of  using  forced  and  proper 
elements  to  describe  their  evolution  in  the  presence  of  drag  force  is  valid.  Figure  3.4 
shows  that  the  proper  inclinations  of  the  dust  particles  are  unchanged  during  their 
evolution  from  1.8  AU  to  0.5  AU. 

In  Figures  3.5  through  3.8  we  compare  the  results  from  the  numerical  integradon 
(shown  as  dots)  and  the  prediction  from  the  dynamical  theory  (solid  lines).  The  forced 
inclination  and  node  will  be  shown  in  Chapter  5  to  be  very  important  parameters 
when  we  construct  the  zodiacal  cloud  model.  After  20,000  years  of  evolution,  the 
difference  in  the  forced  inclination  between  the  dynamical  theory  and  the  actual 
numerical  integration  is  about  0.02  degree.  The  difference  in  the  forced  node  is 
less  than  1  degree.  This  shows  that  the  dynamical  theory  can  describe  properly  the 
evolution  of  dust  particles  under  the  influence  of  both  gravitational  perturbations  and 
PR  drag  as  long  as  the  effects  of  mean  motion  resonance  and  point  mass  scattering 
are  not  important. 

In  Figure  3.9  we  plot  the  forced  and  proper  eccentricity  in  a  log(a)  vs.  log(e) 
frame.  Two  reference  lines  with  slopes  equal  to  1.25  are  also  shown  in  the  diagram. 
As  we  mentioned  in  the  previous  section,  drag  force  is  the  only  factor  that  changes 
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the  proper  eccentricity.  It  changes  the  proper  eccentricity  according  to  Equation  3.19, 
which  predicts  an  1.25  slope  in  the  log(a)  vs.  log(e)  plot.  The  forced  eccentricity 
does  not  approach  the  reference  line  until  the  particles  are  very  close  to  the  Sun.  This 
is  because  both  the  gravitational  perturbation  and  drag  act  on  the  forced  eccentricity. 
When  particles  are  far  away  from  the  Sun,  the  gravitational  perturbation  dominates 
the  change  in  the  forced  eccentricity.  It  is  only  when  the  particles  are  very  far  away 
from  the  planets  (Jupiter  and  Saturn),  that  the  drag  becomes  the  dominant  factor  and 
the  forced  eccentricity  in  the  diagram  approaches  the  straight  line. 

Passage  Through  Resonance 

In  the  real  solar  system,  asteroid  dust  particles  produced  from  the  main  asteroid 
belt  will  go  through  several  Kirkwood  Gaps  and  secular  resonance  regions  while 
they  are  spiralling  towards  the  Sun.  Kirkwood  Gaps  are  the  regions  at  which  a 
particle  is  in  mean  motion  resonance  with  Jupiter  (e.g.  Dermott  and  Murray,  1983). 
For  example,  the  3:1  resonance  is  at  a=2.46  AU  and  the  5:2  resonance  at  2.78 
AU  for  particles  having  /?=0.05037.  The  locations  of  the  Gaps  are  slightly  shifted 
due  to  the  radiation  force.  Secular  resonance  regions  are  the  places  at  which  a 
particle's  orbital  precession  rate  matches  one  or  a  simple  combination  of  several 
of  the  eigenfrequencies  of  the  solar  system.  The  major  secular  resonance  region 
is  located  near  2  AU  (e.g.  Williams,  J.  G.,  1969,  Scholl  et  al.,  1989).  Passage 
through  different  Kirkwood  Gaps  has  some  complicated  effects  on  the  orbits  of  the 
dust  particles.  To  develop  an  analytical  theory  (if  there  is  one)  to  solve  the  problem 
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is  beyond  the  scope  of  this  dissertation.  However,  one  can  observe  their  effects  from 

direct  numerical  integration,  as  we  have  done  here. 

Again,  this  is  a  Sun-Jupiter-Satum  system.  The  initial  forced  elements  are 
calculated  similar  to  the  methods  mentioned  in  the  previous  section.  Two  hundred 
and  fifty  four  Eos  dust  particles  are  placed  at  3.015  AU  at  t=0.  Their  other  orbital 
elements  are  generated  by  the  "particles  in  a  circle"  method  (Dermott  et  al.,  1985, 
1992).  The  changes  in  orbital  elements  due  to  the  releasing  process  are  not  included. 
These  dust  particles  are  integrated  for  50,000  years.  They  end  up  at  about  0.7  AU. 

The  distributions  in  the  (p,q)  and  (h,k)  phase  spaces  of  these  particles  while 
they  are  moving  towards  the  Sun  are  shown  in  Figure  3.10  and  Figure  3.11.  The 
time  interval  between  two  consecutive  panels  is  6,000  years.  The  first  impression 
from  the  figures  is  that  they  remain  in  a  well-defined  circle  in  the  inclination  space 
while  going  through  large  scattering  in  the  eccentricity  space.  In  Figure  3.12  we  plot 
the  proper  inclination  vs.  semimajor  axis  of  the  particles.  The  proper  inclination 
jumps  up  and  down  somewhat  (particularly  at  the  3:1  resonance)  during  the  particles' 
passage  through  Kirkwood  Gaps  even  though  the  magnitude  of  change  is  not  very 
significant  (from  the  observational  point  of  view).  Figure  3.13  and  Figure  3.14 
show  the  variation  in  the  forced  inclination  and  longitude  of  ascending  node  from 
the  numerical  integration  (dots)  and  the  results  from  the  dynamical  theory  (solid 
line).  The  difference  is  quite  small.  This  implies  the  dynamical  theory  works  well 
in  inclination  even  in  the  case  of  mean  motion  resonance  passage.  The  variations 
in  the  forced  eccentricity  and  longitude  of  pericenter  are  in  Figure  3.15  and  Figure 
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3.16.  The  sudden  drop  at  the  3:1  resonance  region  is  very  obvious.  The  difference 
between  the  numerical  results  and  the  dynamical  theory  is  around  30%  near  the  end 
of  the  integration. 

In  all  the  diagrams  there  is  no  dramatic  change  that  can  be  seen  when  particles 
pass  through  the  secular  resonance  region. 

We  have  also  performed  another  test  run  here.  We  changed  the  masses  of  Jupiter 
and  Saturn  to  10%  of  their  real  masses  and  repeated  the  same  numerical  integration. 
This  was  to  reduce  the  strength  of  each  mean  motion  resonance  (Dermott  and  Murray, 
1983).  The  difference  in  the  forced  eccentricity  was  dramatically  reduced  this  time 
as  we  can  see  in  Figure  3.17.  This  also  supports  what  we  have  just  suggested, 
namely  that  the  discrepancies  in  Figures  3.15  and  3.16  are  actually  due  to  passage 
through  mean  motion  resonance  zones. 

From  the  above  results  we  can  conclude  that  passage  through  Kirkwood  Gaps 
does  produce  a  significant  effect  on  particles'  eccentricities  and  longitudes  of  peri- 
center. 

Orbital  Evolution  of  Dust  Particles  from  the  Major  Hirayama  Asteroidal  Families 

The  orbital  evolution  of  dust  particles  in  the  real  solar  system  is  more  complicated 
than  the  cases  presented  previously.  First  of  all,  even  if  the  parent  bodies  which 
produce  dust  particles  constitute  a  well-defined  group  (i.e.  they  form  nice  circles 
in  both  {p,q)  and  {h,k)  spaces),  the  releasing  process  changes  some  of  the  orbital 
elements  of  the  dust  particles  immediately  (see  Equations  3.3  and  3.4).  Thus,  even 
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the  dust  particles  from  a  well-defined  family  do  not  form  a  well-defined  circle  in 
the  eccentricity  space  from  the  very  beginning.  This  makes  the  attempt  to  study 
analytically  the  variation  in  forced  and  proper  eccentricities  of  dust  particles  coming 
from  the  same  family  very  difficult.  Secondly,  inside  the  orbit  of  Jupiter,  there  are 
three  important  planets.  Mars,  Earth,  and  Venus.  When  dust  particles  pass  right  by 
these  planets,  there  are  singularities  (due  to  zero  distance  to  the  planet)  which  cannot 
be  handled  analytically.  The  third  reason  is  that  some  of  the  dust  particles  will  be 
trapped  in  the  near-Earth  mean  motion  resonance  (and  to  a  lesser  extent,  near-Mars 
and  near- Venus  resonances).  Due  to  these  difficulties,  the  practical  way  to  study  the 
orbital  evolution  of  micron-sized  dust  particles  is  direct  numerical  integration. 

In  this  section  we  perform  the  numerical  integration  for  dust  particles  coming 
from  three  major  Hirayama  asteroid  families,  Eos,  Themis,  and  Koronis.  All  the 
planetary  perturbations  (Mercury  and  Pluto  are  not  included),  radiation  pressure, 
Poynting-Robertson  drag,  and  corpuscular  drag  are  included  in  the  calculation.  The 
changes  in  the  orbital  elements  due  to  the  process  of  release  are  also  included. 
Particles  are  assumed  to  leave  their  parent  bodies  with  zero  relative  velocity.  In  order 
to  utihze  the  vector  facility  on  the  IBM  ES/9000,  we  set  up  249  dust  particles  and 
7  planets  for  each  integration.  The  solar  system  is  wound  back  in  time  numerically 
to  an  epoch  such  that  dust  particles  being  released  at  that  epoch  will  reach  a  chosen 
heliocentric  distance  in  1983.  The  initial  forced  elements  of  dust  particles  are 
calculated  from  the  orbital  elements  of  the  7  planets  at  the  beginning  epoch.  The 
proper  nodes  and  longitudes  of  pericenter  of  dust  particles  are  chosen  randomly 
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between  0°  and  360°.  The  proper  eccentricities  and  inclinations  of  families  are  from 

Dermott  et  al.    (1985). 

9  /fm  Diameter  Particles 

We  first  integrate  one  wave  of  particles  from  each  family.  Each  wave  ends  up 
inside  1  AU  in  1983.  For  Eos  particles,  the  variations  in  (p,q)  and  {h,k)  spaces  during 
their  evolution  are  shown  in  Figures  3.18  and  3.19.  The  first  panel  in  each  figure 
is  the  distribution  of  the  parent  bodies  which  produce  the  dust  particles.  The  time 
interval  between  two  consecutive  panels  is  6,500  years.  The  label  above  each  panel 
is  the  average  semimajor  axis  of  the  dust  particles  at  that  time. 

In  the  inclination  space,  Eos  dust  particles  remain  in  a  circle  even  after  they  pass 
the  Earth  while  Themis  and  Koronis  dust  particles  have  large  scatterings  when  they 
approach  the  Earth.  This  is  because  Themis  and  Koronis  dust  particles  have  much 
smaller  inclinations  compared  to  that  of  Eos  dust  particles.  The  probability  of  low 
inclination  objects  being  scattered  by  terrestrial  planets  is  higher  than  that  of  high 
inclination  objects.  In  the  eccentricity  space,  Poynting-Robertson  drag  reduces  the 
eccentricities  of  dust  particles  while  they  are  moving  towards  the  Sun.  However, 
when  they  approach  the  Earth  some  of  them  are  trapped  in  the  near-Earth  mean 
motion  resonances.  These  particles  gain  energy  and  lose  angular  momentum  during 
the  process.  Eventually  they  will  escape  from  the  resonance  zone  and  continue  their 
journey  towards  the  Sun.  This  phenomenon  will  be  discussed  in  more  detail  in  the 
next  section. 
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Obtaining  the  particle  distribution  at  different  heliocentric  distances  in  1983  (for 
the  purpose  of  constructing  the  zodiacal  cloud  model)  is  different  from  studying  the 
dynamical  evolution  of  the  particles.  Different  waves  (each  "wave"  is  the  numerical 
integration  of  the  orbits  of  249  dust  particles)  of  particles  must  be  sent  out  at  different 
epochs  such  that  they  end  up  at  different  heliocentric  distances  in  1983.  This  task 
requires  considerable  CPU  time.  The  total  data  for  three  families  took  about  3 
months  on  an  IBM  ES/9000  supercomputer.  For  Eos  dust  particles,  we  sent  in  23 
different  waves  dating  back  to  the  year  57,500  years  before  1983.  The  time  interval 
between  each  wave  is  2,500  years.  Then  we  analyzed  the  data  at  the  end  of  each 
integration  and  obtained  the  orbital  elements  distribution  at  different  heliocentric 
distances.  Thus,  we  have  the  forced  elements  as  a  function  of  heliocentric  distance 
from  around  0.5  AU  to  3.2  AU.  The  forced  inclinations  and  nodes  for  Eos  particles 
are  shown  in  Figures  3.20  and  3.21.  Themis  and  Koronis  particles  have  similar 
forced  inclinations  and  nodes  until  they  reach  around  1  AU  where  they  no  longer 
remain  in  a  circle  due  the  gravitational  scattering  by  the  Earth.  In  the  (h,k)  space, 
particles  are  scattered  even  from  the  very  beginning  as  we  can  see  in  Figure  3.19. 
Hence,  we  did  not  attempt  to  find  the  corresponding  forced  elements  by  fitting  a 
circle  to  the  points  in  eccentricity  space.  Instead,  we  used  a  tabular  form  from  the 
data  to  generate  the  eccentricities  and  longitudes  of  pericenter  for  the  dust  particles 
when  we  constructed  the  zodiacal  cloud  model. 

Once  the  orbital  element  distributions  of  dust  particles  from  a  known  family  have 
been  obtained,  it  is  very  easy  to  construct  a  real  three  dimensional  model  zodiacal 
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cloud  from  the  results.  We  did  just  that  for  Eos  dust  particles.  Figure  3.22  is  the 
vertical  cross-section  picture  of  the  Eos  cloud.  This  is  the  cloud  which  may  be 
responsible  for  the  10°  dust  band.  Actually  the  so-called  "solar  system  dust  band" 
is  not  a  band.  Because  of  the  drag  effects,  dust  particles  go  all  the  way  from  their 
origin  to  the  Sun  (the  reason  there  is  a  hole  in  the  figure  inside  0.5  AU  is  that  we 
terminate  the  integration  when  dust  particles  reach  there).  The  actual  "band"  is  the 
enhancement  of  the  observed  flux  at  certain  ecliptic  latitudes  due  to  the  concentrations 
of  dust  particles  at  the  edges.  This  concentration  can  be  seen  in  Figure  3.22.  In  the 
real  solar  system,  this  cloud  is  embedded  in  the  zodiacal  cloud  which  includes  dust 
particles  from  many  different  sources.  However,  this  concentration  is  strong  enough 
such  that  when  the  whole  zodiacal  cloud  was  scanned  by  IRAS,  this  "band"  (more 
precisely,  a  tiny  peak  on  top  of  the  background)  at  around  10°  was  observed. 

The  forced  inclination  and  longitude  of  the  ascending  node  of  the  dust  particles 
determine  the  plane  of  symmetry  of  the  zodiacal  cloud.  This  plane  is  a  warped  plane 
as  the  forced  inclination  and  longitude  of  the  ascending  node  vary  with  heliocentric 
distances.  When  we  observe  the  zodiacal  cloud,  the  latitudes  of  the  peak  fluxes 
around  the  sky  are  determined  by  this  warped  plane.  These  positions  are  very 
well-defined  in  the  IRAS  observations.  One  of  the  most  important  achievements 
of  the  present  work  is  being  able  to  use  the  forced  inclination  and  longitude  of 
the  ascending  node  from  the  dynamical  study  and  the  SIMUL  model  to  predict  the 
latitudes  of  the  peak  fluxes  as  observed  by  IRAS.  Comparison  of  those  predictions 
with  the  observations  will  be  presented  in  Chapter  5. 


41 
In  Figure  3.20  the  data  show  that  the  forced  inclination  increases  towards  small 
heliocentric  distances.  This  increase  had  been  reported  from  observations  prior  to 
the  launch  of  IRAS  (e.g.  Misconi  and  Weinberg,  1978,  Leinert  et  al.,  1980).  Is  this 
increase  due  to  terrestrial  planets,  especially  Venus  (which  has  a  very  high  inclination, 
3.4°)  as  some  early  papers  suggested  (Misconi  and  Weinberg,  1978,  Gustafson  and 
Misconi,  1986)?  We  examine  this  by  taking  out  the  terrestrial  planets  one  at  a  time 
and  integrating  the  orbits  of  Eos  dust  particles  in  each  case.  The  resulting  forced 
inclinations  from  these  cases  are  in  Figure  3.23.  The  effect  of  the  terrestrial  planets 
can  be  seen  clearly.  The  increase  in  the  forced  inclination  has  nothing  to  do  with  the 
existence  of  the  terrestrial  planets.  It  is  the  passage  through  the  secular  resonance 
zone  that  increases  the  forced  inclination. 

4  fxm  and  14  //m  Diameter  Particles 

We  assume  that  9  fim  diameter  dust  particles  are  the  major  contributor  to  the  25 
//m  waveband  observation  of  IRAS.  So  far  we  have  studied  the  dynamical  evolution 
of  9  fxm  diameter  dust  particles  and  their  distribution  in  the  solar  system.  In  Chapter 

5  we  will  show  a  zodiacal  cloud  model  based  on  particles  of  this  one  size.  However, 
particles  of  different  sizes  do  contribute  to  the  25  /im  waveband  observations  as 
well.  Eventually  the  model  zodiacal  cloud  must  include  results  from  particles  with 
a  range  of  sizes.  In  this  section  we  present  the  dynamical  study  for  particles  with 
diameters  equal  to  4  /im  and  14  /xm.  This  covers  a  range  of  3.5  in  sizes  (a  factor 
of  43  in  masses  if  they  have  the  same  density). 
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The  integration  procedures  here  are  similar  to  those  in  the  previous  section.  The 
corresponding  CPU  time  for  14  /^m  diameter  dust  particles  of  course  is  much  longer 
than  that  of  9  /im  diameter  dust  particles. 

The  forced  inclinations  and  nodes  as  functions  of  semimajor  axes  for  the  three 
different  size  particles  are  shown  in  Figures  3.24  and  3.25.  Among  them,  the  14  ;im 
diameter  dust  particles  have  the  smallest  drag  rate.  They  pick  up  the  most  increase 
in  the  forced  inclination  while  they  pass  through  secular  resonance  regions  and  end 
up  with  the  largest  forced  inclination  in  the  inner  solar  system.  The  difference  in 
inclination  at  1  AU  between  these  three  different  size  particles  is  about  1.5  degree. 
Obviously  this  difference  must  be  considered  when  a  better  zodiacal  cloud  model 
(comprising  particles  of  4,  9,  and  14  /im  in  diameters)  is  to  be  constructed.  How  to 
combine  these  results  is  very  important  in  constructing  a  successful  zodiacal  cloud 
model  in  the  future. 

Near-Earth  Dust  Grains 

In  this  section  we  will  describe  some  properties  of  asteroidal  dust  particles  when 
they  are  in  the  vicinity  of  the  Earth.  These  include  seasonal  variations  in  the  number 
of  particles  that  intersect  the  Earth,  how  to  identify  the  origins  of  near-Earth  dust 
particles,  and  the  phenomenon  of  trapping  in  near-Earth  mean  motion  resonances. 

The  forced  inclination  and  longitude  of  ascending  node  determine  the  plane  of 
symmetry  of  the  zodiacal  cloud.  They  also  produce  another  interesting  phenomenon 
for  particles  with  small  proper  inclinations:  the  seasonal  variation  of  dust  particles 
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collected  in  near-Earth  orbit.  The  distributions  in  inclination  space  for  4  /im  and 
9  i^im  dust  particles  are  shown  in  Figure  3.26.  These  are  particles  in  the  Earth- 
crossing  region  in  1983  found  from  our  numerical  integrations.  Notice  the  off-center 
distribution  for  Themis  and  Koronis  particles.  There  are  more  Themis  and  Koronis 
particles  in  the  first  and  fourth  quadrants  than  in  the  second  and  third  quadrants  of 
the  diagram.  When  the  Earth  moves  around  its  orbit  over  one  year,  it  will  encounter 
more  Themis  and  Koronis  dust  particles  at  certain  places.  These  places  correspond 
to  the  areas  in  the  inclination  space  where  there  are  more  particles  (first  and  fourth 
quadrants).  This  kind  of  seasonal  variation  for  4  and  9  /<m  particles  is  shown  in 
Figure  3.27.  Here  we  have  a  histogram  showing  the  variation  in  the  relative  number 
of  particles  as  a  function  of  the  time  of  the  year.  The  Earth  encounters  more  Themis 
particles  at  their  ascending  nodes  in  November  than  in  April.  The  difference  is  more 
than  a  factor  of  8. 

There  is  no  significant  seasonal  variation  for  particles  with  large  proper  inclina- 
tions, like  Eos.  In  the  vicinity  of  the  Earth  the  proper  inclination  of  Eos  particles 
is  much  larger  than  the  forced  inclination.  Consequently,  even  though  the  circle  is 
shifted  from  the  origin  in  the  inclination  space  (see  Figure  3.26),  the  difference  in 
terms  of  number  of  particles  in  each  quadrant  is  not  very  large.  Hence  the  seasonal 
variation  for  Eos  particles  in  Figure  3.27  is  not  as  dramatic  as  that  of  Themis  particles. 

Themis,  Koronis,  and  Eos  are  the  most  prominent  families  in  the  asteroid  belt. 
The  ratio  of  dust  particles  produced  from  these  three  families  to  non-family  asteroids 
is  1  :    2.5.    The  ratio  of  dust  particles  produced  among  these  three  families  is 
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Themis  :  Eos  :  Koronis  =  1  :  0.61  :  0.12  (Durda,  1993).  This  means  Themis 
and  Eos  produce  17%  and  10%  of  all  the  asteroidal  dust  particles,  respectively. 
This,  together  with  the  facts  that  (1)  particles  from  certain  families  do  not  arrive 
uniformly  in  time  (e.g.  most  Themis  particles  with  ascending  nodes  encounter  the 
Earth  around  November)  and  (2)  particles  from  one  family  have  a  well-defined  range 
of  their  inclinations,  makes  it  possible  to  identify  the  origins  of  some  family  dust 
particles.  In  Figure  3.28  we  show  the  relative  numbers  of  particles  as  a  function  of 
inclination.  These  are  4  and  9  ^m  dust  particles  at  their  ascending  nodes  which  will 
encounter  the  Earth  in  November.  Themis  and  Eos  particles  are  represented  by  bold 
lines  in  the  histogram  while  the  whole  non-family  asteroidal  particles  are  shown  as 
thin  lines.  Much  useful  information  can  be  obtained  from  this  figure.  First  of  all, 
the  low  inclination  dust  particles  in  November  are  dominated  by  Themis  members. 
Half  of  the  particles  around  12°  inclination  are  of  Eos  origin.  If  we  utilize  these 
results  to  interpret  the  samples  from  an  Earth-orbiting  dust-collection  probe,  we  can 
obtain  detailed  information  (e.g.  chemical  composition)  about  some  asteroid  families 
without  having  to  go  to  the  asteroid  belt.  For  example,  from  this  dynamical  study  we 
know  that  half  of  the  dust  particles  being  collected  in  November  with  12°  inclination 
are  Eos  particles.  Then,  we  can  collect,  say,  20  particles  with  inclinations  around 
12°  and  analyze  their  composition.  Half  of  them  will  have  similar  compositions 
(if  we  assume  particles  from  the  same  family  have  the  same  compositions).  These 
must  be  Eos  particles.  Their  chemical  compositions  are  those  of  the  Eos  asteroid 
family.  Also,  most  of  the  low  inclination  (around  3°)  dust  particles  being  collected 
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in  November  must  have  similar  compositions.  These  are  particles  from  the  Themis 

family. 

An  important  objective  for  any  future  dust  collection  mission  is  to  establish  a 
relationship  between  the  dust  particles  being  collected  and  their  parent  bodies.  Our 
results  here  have  shown  that  it  may  be  possible  to  achieve  such  a  goal.  With  modem 
techniques,  intact  capture  of  particles  with  velocities  less  than  9  km/sec  is  possible 
(these  are  most  likely  to  be  asteroidal  dust  particles).  From  our  study,  we  can  provide 
not  only  the  best  time  but  also  the  orientation  of  the  orbits  of  the  dust  particles  that 
can  be  collected  and  related  to  several  known  asteroidal  families. 

When  dust  particles  approach  the  Earth  from  outside,  some  of  them  will  be 
trapped  in  the  near-Earth  mean  motion  resonances.  The  trapping  probability  depends 
on  many  factors.  These  include  the  drag  rate,  the  orbital  elements  of  the  particles, 
the  strength  of  the  resonance,  and  the  geometry  of  how  a  particle  approaches  the 
resonance.  This  trapping  in  mean  motion  resonance  will  produce  some  interesting 
observational  consequences  (e.g.  Jackson  and  Zook,  1989,  Reach,  1991,  Dermott  et 
al.,  1993).  Here  we  will  describe  qualitatively  the  behavior  of  particles  when  they 
approach  the  Earth  based  on  our  numerical  results.  Similar  phenomena  happen  near 
Venus  and  Mars,  but  to  a  lesser  extent. 

The  orbital  evolution  of  particles  being  trapped  in  a  mean  motion  resonance  with 
a  planet  is,  actually,  very  easy  to  understand.  When  a  particle  is  spiralling  towards  the 
Sun,  the  drag  force  takes  away  its  energy  and  angular  momentum.  When  a  particle 
is  trapped  in  a  resonance,  it  gains  energy  from  the  resonance  to  counterbalance  the 
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loss  due  to  drag  force.  Whether  a  particle  will  be  trapped  or  not  depends  on  how 
much  energy  it  gains  when  it  passes  through  a  given  resonance  zone.  If  it  does 
not  obtain  enough  energy,  it  simply  passes  right  through  (with  some  minor  changes 
in  its  orbital  elements).  Otherwise  it  will  get  trapped.  However,  this  particle  will 
not  be  trapped  forever.  Eventually  it  will  have  a  series  of  close  encounters  or  a 
very  deep  close  encounter  with  the  planet.  Then,  the  close  encounter  removes  this 
particle  from  the  resonance. 

In  our  numerical  integration,  we  have  studied  the  evolution  of  particles  with 
diameters  equal  to  4,  9,  and  14  /xm.  The  drag  rate  for  small  dust  particles  is  larger 
than  that  for  the  big  particles.  Hardly  any  of  the  4  //m  dust  particles  are  trapped  in 
any  resonances.  A  few  of  the  9  ^m  dust  particles  are  trapped.  In  Figure  3.29  we 
show  a  wave  of  9  fim  Eos  particles  moving  towards  the  Earth.  In  Figure  3.30  we 
show  the  semimajor  axes  vs.  time  for  50  of  them  when  they  approach  the  Earth. 
Some  of  them  do  get  trapped,  but  the  overall  trapping  in  resonance  is  not  very 
significant.  For  14  /im  dust  particles,  the  trapping  percentage  is  higher  because  of 
their  lower  drag  rate. 
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CHAPTER  4 
ORBITAL  EVOLUTION  OF  COMETARY  PARTICLES 

The  major  difficulty  in  studying  the  dynamical  evolution  of  cometary  dust 
particles  and  their  contribution  to  the  zodiacal  cloud  is  that  of  identifying  their 
sources.  Unlike  asteroidal  dust  particles,  which  are  mostly  from  the  main  asteroid 
belt,  cometary  dust  particles  have  no  well-defined  parent  bodies.  We  do  not  know 
whether  they  are  from  a  highly  active  comet,  a  group  of  comets,  or  all  short  period 
comets.  This  makes  the  attempt  to  find  the  contribution  of  cometary  particles  to  the 
zodiacal  cloud  seem  impossible.  Nevertheless,  we  feel  that,  even  if  we  don't  know 
the  exact  source  of  cometary  dust  particles,  by  studying  the  characteristics  of  the 
model  zodiacal  cloud  produced  from  a  typical  cometary  source,  we  will  gain  some 
insight  into  the  overall  problem. 

The  typical  cometary  source  we  choose  is  comet  Encke.  In  this  chapter  we 
study  the  orbital  evolution  of  Encke-type  dust  particles  numerically  and  find  their 
distribution  in  the  solar  system.  We  also  show  that  comets  cannot  be  the  source  of 
the  solar  system  dust  bands  observed  by  IRAS. 

Orbital  Evolution  of  Encke-type  Dust  Particles 

Encke  has  long  been  proposed  to  be  the  major  source  of  the  zodiacal  cloud 
(e.g.    Whipple,  1967).    Encke  is  special  in  the  sense  that  it  is  dynamically  stable 
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against  the  gravitational  perturbation  of  Jupiter  (Opik,  1963)  and  because  of  its  strong 
correlation  with  several  meter  showers  and  the  so-called  Taurid  Complex  (e.g.  Steel, 
et  al,  1991a,b,  Hartung,  1993).  It  is  thus  natural  to  study  the  orbital  evolution  of 
Encke-type  dust  particles  when  there  are  no  other  better  choices.  Our  definition 
for  Encke-type  dust  particles  is  those  with  initial  semimajor  axes,  eccentricities  and 
inclinations  similar  to  those  of  Encke.  The  mean  semimajor  axis,  eccentricity,  and 
inclination  for  Encke  are  2.2  AU,  0.85,  and  12°,  respectively.  The  initial  longitudes 
of  ascending  nodes,  longitudes  of  pericenters,  and  mean  longitudes  of  dust  particles 
are  randomly  chosen  between  0°  and  360°.  This  means  these  micron-sized  dust 
particles  are  not  released  by  means  of  direct  sublimation  from  the  parent  body  when 
it  is  near  its  pericenter.  The  PR  drag  lifetime  for  Encke-type  micron-sized  dust 
particles  is  about  5,000  years.  If  they  are  produced  from  a  comet  at  its  pericenter, 
they  will  not  have  time  to  precess  and  spread  out  their  orbits  to  produce  the  axial 
symmetry  structure  of  the  zodiacal  cloud.  We  assume  that  it  is  the  disruption  of 
large  Encke-type  cometary  meteoroids  (large  enough  such  that  their  longitudes  of 
ascending  nodes  and  longitudes  of  pericenters  have  been  dispersed  between  0°  and 
360°)  that  produces  these  dust  particles.  Griin  et  al.  (1985a,b)  concluded  that  this  is 
the  most  efficient  way  to  deposit  cometary  dust  particles  into  the  zodiacal  cloud. 

We  first  set  up  the  initial  condition  for  five  hundred  9  i.im  diameter  dust  particles 
by  the  "particles  in  a  circle"  method  (Dermott  et  al.,  1985,  1992).  We  start  the 
integration  at  five  thousand  years  before  1983.  Dust  particles  are  then  integrated  for 
five  thousand  years.   We  record  their  orbital  elements  once  every  100  years.   The 


79 
variations  in  the  (h,k)  and  (p,q)  spaces  when  the  particles  are  spiralling  towards  the 
Sun  are  shown  in  Figures  4.1  and  4.2.  The  time  difference  between  two  consecutive 
panels  is  six  hundred  years.  The  first  panel  is  the  distribution  of  the  parent  bodies 
at  time  equal  to  zero. 

The  major  difference  between  Figure  4.2  and  Figure  3.18  is  that  cometary  dust 
particles  do  not  remain  in  a  well-defined  circle  in  inclination  space  right  after  they  are 
released.  If  we  examine  the  inclination  of  each  individual  particle,  we  will  find  that 
it  goes  through  huge  variations  as  the  particle  moves  towards  the  Sun.  As  a  whole, 
cometary  particles  form  a  doughnut-like  distribution  with  well-defined  maxima  and 
minima  in  the  inclination  space.  These  features  are  the  most  important  characteristics 
of  cometary  type  dust  particles.  They  may  play  crucial  roles  in  combining  with  the 
asteroidal-type  dust  particles  to  form  the  zodiacal  cloud.  This  will  be  discussed  in 
detail  in  Chapter  5.  The  physical  interpretation  for  these  inclination  variations  is 
discussed  in  the  next  section. 

The  evolution  of  Encke-type  dust  particles  in  (a,e)  space  is  shown  in  Figure 
4.3.  These  five  groups  correspond  to  data  from  one  thousand  (large  filled  circles) 
to  five  thousand  (crosses)  years  after  they  were  released.  Very  few  particles  are 
trapped  in  mean  motion  resonances  with  Jupiter  or  the  Earth.  This  is  due  to  their 
high  eccentricities  which  give  them  much  higher  drag  rates  than  asteroidal  particles 
with  the  same  /?  and  thus  very  much  lower  probabilities  of  capture.  We  demonstrate, 
in  Figure  4.4,  the  relationship  between  the  drag  rate,  da/dt,  and  the  eccentricity  of  a 
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particle.  The  drag  rate  of  the  dust  particle  is  given  by  (see  Equation  3.13) 

da  const.{2 -\- "ie^)  const..,  ,  ^^  ,^ 

:j7  = ^ r^  = /(e)  (4.1) 

«^  a(l  -  e2)2  « 

where  a  and  e  are  the  semimajor  axis  and  eccentricity  of  the  dust  particle,  respec- 
tively. With  the  same  starting  semimajor  axis,  the  drag  rate  of  the  dust  particle 
increases  dramatically  towards  the  high  eccentricity  range.  For  example,  da/dt  of  an 
eccentricity  0.85  particle  is  nearly  14  times  larger  than  that  of  an  eccentricity  0.1 
particle.  Comparing  Encke-type  dust  particles  with  typical  asteroidal  dust  particles, 
we  notice  that  this  high  drag  rate  makes  them  move  much  faster  through  the  mean 
motion  resonance  zone  in  such  a  way  that  they  do  not  have  time  to  gain  enough 
energy  from  the  resonance  to  counterbalance  the  loss  due  to  PR  drag  orbital  decay. 

Gauss'  Equations 

The  reason  why  Encke-type  dust  particles  go  through  huge  variations  in  their 
inclinations  can  be  understood  from  Gauss'  equations.  These  equations  provide  a 
better  intuitive  understanding  of  the  gravitation  perturbations  from  planets  (see,  for 
example,  Danby,  1988)  than  the  Lagrange  equations.  The  rate  of  change  of  the 
inclination  of  the  particle  is  given  by 

dl  r 

—  =  ,W  cos  {io  +  /)  (4.2) 

dt        na-\/\  —  e- 

where  r  is  the  heliocentric  distance  of  the  dust  particle,  to  and  /  are  the  argument 
of  pericenter  and  the  true  anomaly  of  the  dust  particle,  respectively.  W  is  the 
planetary  perturbation  force  component,  acting  perpendicularly  to  the  orbital  plane 
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of  the  particle.  When  we  consider  the  long-term  average  perturbation,  it  is  obvious 
that  as  the  eccentricity  of  a  dust  particle's  orbit  increases,  both  the  factor  -7=^  and 
W  increase.  This  will  lead  to  an  increase  not  only  in  the  rate  of  change  of  inclination 
but  also  in  the  maximum  and  minimum  inclinations  that  particle  can  reach.  Particles 
starting  with  the  same  orbital  elements  will  reach  the  same  maxima  and  minima  in 
their  inclinations. 


Can  Comets  be  the  Source  of  the  Solar  System  Dust  Bands? 

IRAS  discovered  several  dust  bands  during  its  all-sky  survey  in  1983.  The 
central  band  (later  it  was  found  that  there  are  actually  two  central  bands)  and  the  10° 
band  are  very  clear  in  all  four  wavelengths.  Several  theories  have  been  proposed 
to  explain  their  origins.  Dermott  et  al.  (1984,  1985)  first  pointed  out  the  possible 
correlation  between  these  bands  and  three  major  Hirayama  asteroidal  families  and 
later  confirmed  it  from  their  dynamical  study  and  modeling  results  (Dermott  et  al., 
1990,  1992,  1993).  The  central  bands  are  probably  due  to  the  Themis  and  Koronis 
families  while  the  Eos  family  may  be  responsible  for  the  10°  band.  However,  the 
possibility  still  exists  that  some  comets  may  produce  these  bands  (Dermott  et  al., 
1984,  Clube  and  Asher,  1990).  Among  all  the  possible  candidates,  Encke  has  been 
proposed  to  be  the  source  of  the  10°  band.  This  is  primarily  due  to  the  fact  that 
Encke 's  mean  inclination  is  very  close  to  10  degrees.  Here  we  will  discuss,  from 
the  dynamical  point  of  view,  why  a  comet  cannot  be  the  source  of  any  dust  bands. 
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Before  we  show  why  there  is  no  cometary  origin  of  dust  bands,  we  must  first 
understand  how  to  form  a  dust  band. 

In  the  (p,q)  phase  space,  the  osculating  inclination  of  a  dust  particle  is  the 
vectorial  sum  of  two  components,  the  forced  and  the  proper  components.  The 
proper  inclinations  of  asteroidal  dust  particles  remain  unchanged  during  their  orbital 
evolutions  (see  Figure  3.18,  this  is  a  key  feature  to  form  a  dust  band)  while 
their  proper  nodes  precess  about  a  point  in  inclination  space  defined  by  the  forced 
inclination  and  node.  From  the  geometric  point  of  view,  this  means  that  the  orbits 
of  asteroidal  dust  particles  precess  about  the  axis  of  a  common  plane,  the  plane 
of  symmetry,  which  is  defined  by  the  forced  inclination  and  node.  Dust  particles 
produced  from  the  same  family  will  form  a  doughnut-like  structure  if  the  family 
members  are  old  enough  such  that  their  nodes  have  been  spread  all  over  the  sky. 
This  doughnut-like  structure  is  tilted  to  the  ecliptic  with  the  inclination  and  node 
defined  by  the  forced  inclination  and  node.  When  we  observe  this  structure  from 
the  Earth,  we  will  see  two  maxima  in  flux  at  the  two  edges  of  the  doughnut  at  any 
given  ecliptic  longitude.  This  is  the  so-called  "dust  band".  These  maxima  still  exist 
even  when  we  consider  the  fact  that  dust  particles  are  spiralling  towards  the  Sun  and 
fill  up  the  space  between  their  origin  and  the  Sun  (see  Figure  3.22).  With  this  band- 
forming  mechanism  in  mind,  we  can  now  examine  the  orbital  evolution  of  cometary 
dust  particles  and  see  if  they  are  capable  of  forming  a  dust  band. 

As  we  mentioned  previously,  while  asteroidal  dust  particles  remain  in  a  well- 
defined  circle  in  the  inclination  space,  cometary  dust  particles  go  through  huge 
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variations  in  their  inclinations.  From  the  geometric  point  of  view,  this  means  there 

will  be  no  well-defined  doughnut  like  structure,  hence  no  well-defined  edges.  Thus, 
it  is  impossible  for  cometary  particles  to  form  a  dust  band. 

Another  observational  consequence  from  this  huge  variation  in  inclination  is  that 
cometary  dust  particles  will  contribute  more  (relatively  speaking,  when  compared 
with  asteroidal  dust  particles)  to  the  observed  flux  at  high  ecliptic  latitudes.  In 
Chapter  5,  we  will  show  that  this  may  be  the  key  feature  needed  to  account  for  the 
observed  shape  of  the  zodiacal  cloud. 

Distribution  of  Encke-type  Dust  Particles  in  the  Solar  System 

In  order  to  find  the  contribution  of  Encke-type  dust  particles  to  the  zodiacal 
cloud  we  must  have  an  idea  of  the  distribution  of  these  particles  in  the  solar  system 
at  the  time  of  the  IRAS  observations.  We  need  the  distributions  of  all  the  orbital 
elements.  We  have  constructed  two  models,  to  be  described  in  this  section.  The  only 
difference  between  these  two  models  is  in  their  initial  inclinations.  The  procedures 
used  to  obtain  the  dust  distributions  are  similar  to  those  in  Chapter  3.  We  send 
out  waves  of  particles  once  every  300  years  starting  from  the  year  5,400  years 
before  1983.  Most  of  the  particles  being  released  before  the  year  (1983-5,400) 
will  have  semimajor  axes  less  than  0.1  AU  by  the  year  1983  and  thus  make  no 
contributions  to  the  IRAS  flux  data.  Each  wave  contains  500  dust  particles.  Their 
orbits  are  integrated  numerically.  All  the  planetary  perturbations  from  7  planets, 
radiation  pressure,  Poynting-Robertson  drag,  and  corpuscular  drag  are  included  in 
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the  numerical  integrations.  The  (3  for  the  dust  particles  is  assumed  to  be  0.05037. 
The  corpuscular  drag  is  assumed  to  be  35  %  that  of  PR  drag.  The  initial  semimajor 
axes  and  proper  eccentricities  of  particles  are  2.2  AU  and  0.85,  respectively.  In  the 
first  model,  the  initial  proper  inclinations  of  all  the  dust  particles  are  set  equal  to 
Encke's  mean  inclination,  12°.  In  the  second  model,  the  initial  proper  inclinations 
of  the  particles  are  generated  according  to  Encke's  variation  in  inclination  in  the 
past.  The  longitudes  of  ascending  node,  the  longitudes  of  pericenter,  and  the  mean 
longitudes  of  the  dust  particles  are  generated  randomly  between  0°  and  360°.  The 
forced  elements  at  2.2  AU  are  calculated  according  to  the  orbital  elements  of  the 
planets  at  the  time  of  particle  release.  All  the  calculations  are  performed  on  an  IBM 
ES/9000.  We  exclude  particles  when  their  semimajor  axes  become  less  than  0. 1  AU 
or  when  their  orbits  become  unbounded.  For  each  integration,  all  six  orbital  elements 
(a,  e,  I,  n,  €c,  A)  of  the  dust  particles  are  recorded  every  100  years  until  the  end  of 
the  integration.  Finally  we  combine  the  results  from  the  end  of  each  integration  to 
get  the  orbital  elements  of  particles  at  different  heliocentric  distances  in  1983. 

The  final  distributions  of  dust  particles  from  the  second  model  (nearly  8,000 
of  them)  in  (a,e),  (a,I),  and  (Q,u7)  phase  spaces  are  shown  in  Figures  4.5  to  4.7. 
Few  particles  are  found  to  be  trapped  in  mean  motion  resonances  with  Jupiter  or 
the  Earth,  as  we  expected. 

Figure  4.6  shows  the  distribution  in  inclination  as  a  function  of  semimajor  axis 
in  1983.  Particles  are  distributed  quite  uniformly  between  2°  and  35°.  Particles  with 
high  inclinations  will  contribute  more  to  the  high  ecliptic  latitude  flux. 
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The  relationship  between  U  and  uj  of  cometary  dust  particles  is  an  interesting 
one,  as  we  can  see  in  Figure  4.7.  Most  of  the  particles  seem  to  have  their  arguments 
of  pericenters,  lj,  concentrated  near  0°  and  180°.  To  understand  this,  we  need  to  use 
Gauss'  equations  again.  It  can  be  easily  shown  (e.g.  see  Whipple,  1940)  that  the  rate 
of  change  of  w  of  a  high  eccentricity  particle  under  the  gravitational  perturbation  of 
Jupiter  is 

^  =  ^^1-^-2  cos  (2c^)  (4.3) 

where  kj  and  k2  are  positive  constants. 

From  the  above  equation,  we  can  see  that  ^  has  minima  at  a;  equal  to  0°  and 
180°,  and  has  maxima  at  u  equal  to  90°  and  270°.  Thus  the  orbit  of  a  cometary  dust 
particle  spends  most  of  its  time  around  u  equal  to  0°  and  180°.  That  is  the  reason 
why  there  are  diagonal  structures  in  Figure  4.7. 

We  next  divide  all  the  particles  into  groups  according  to  the  different  ranges  of 
their  semimajor  axes  and  plot  their  distributions  in  the  (h,k)  and  (p,q)  spaces  in  Figure 
4.8  and  Figure  4.9.  Each  panel  in  the  figures  contains  particles  being  released  from 
different  epochs.  However,  in  the  (p,q)  space  the  overall  distribution  of  particles  is 
quite  similar;  they  spread  out  in  a  disk-like  structure  with  maximum  radius  roughly 
equal  to  35°.  In  the  (h,k)  space,  the  particles  do  not  behave  like  asteroidal  dust 
particles  (see  Figure  3.19,  where  they  are  scattered  all  over  the  phase  space).  There 
are  several  reasons  for  this.  First  of  all,  the  changes  in  the  eccentricities  when 
dust  particles  are  released  from  their  parent  bodies  is  small  compared  with  the  high 
eccentricity  of  the  parent  comet,  0.85.  Secondly,  because  of  their  large  eccentricities, 
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mean  motion  resonances  (that  is,  passage  through  mean  motion  resonance  and 
trapping  in  mean  motion  resonance)  do  not  affect  their  eccentricities  significantly. 
Also,  because  the  inclinations  of  cometary  dust  particles  go  through  huge  variations 
(in  this  case,  from  2°  to  35°),  the  planetary  scattering  effect  is  also  insignificant. 

The  distributions  of  the  orbital  elements  of  particles  from  the  first  model  are  very 
similar  to  the  results  from  the  second  model.  The  only  difference  is  in  their  inclination 
distribution.  Figure  4.10  shows  the  ip,q)  distributions  at  different  semimajor  axis 
ranges  from  the  first  model.  If  we  compare  this  figure  with  Figure  4.9,  we  can  see 
that  in  this  diagram  there  is  an  empty  hole  near  the  origin  in  each  panel.  Also,  the 
maximum  inclination  of  the  particles  is  only  around  30°.  As  we  mentioned  earlier, 
the  first  model  starts  the  particles  with  12°  proper  inclinations  while  the  second  model 
generates  the  proper  inclinations  of  dust  particles  by  the  variation  in  inclination  of 
Encke  in  the  past.  Our  backward  numerical  integration  of  Encke  shows  that  it  spent 
more  time  around  5°  and  15°  than  near  10°.  Thus  there  are  more  dust  particles 
in  the  second  model  having  both  smaller  and  higher  (than  12°)  proper  inclinations 
when  the  integration  begins.  It  is  the  variation  in  inclination  of  the  small  inclination 
particles  (initially  with  proper  inclinations  around  5°)  that  fill  up  the  hole  in  the  {p,q) 
space,  as  shown  in  Figure  4.10.  Particles  starting  with  higher  inclinations  (around 
15°)  push  the  maximum  inclination  to  around  35°.  The  difference  in  their  inclination 
distribution  produces  significant  differences  in  the  model  zodiacal  clouds.  Figures 
4.1 1  and  4.12  show  the  distributions  of  inclinations  in  1983  from  the  first  and  second 
model,  respectively.  There  is  an  increase  in  terms  of  relative  numbers  towards  the 
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peak  at  about  25°  from  the  first  model,  compared  to  the  peak  at  around  8°  from 
the  second  model.   From  the  observational  point  of  view  this  means  that  the  first 
model  produces  a  zodiacal  cloud  background  with  much  higher  wings  at  high  ecliptic 
latitudes  than  the  second  model.  We  will  discuss  this  in  more  detail  in  Chapter  5. 

The  program  which  calculates  the  observable  flux  from  the  particle  distribution 
is  called  SIMUL  (Dermott  et  al.,  1992,  Xu  et  al.,  1993).  The  SIMUL  program 
first  generates  orbits  at  different  semimajor  axes  according  to  the  number  density 
function,  n(a),  then  generates  the  other  four  orbital  elements  of  the  orbits,  e,  I,  0, 
and  w  at  any  given  semimajor  axis  using  the  results  from  the  dynamical  study.  The 
number  of  orbits  per  unit  volume  at  a  given  a  is  defined  by  n(a).  For  asteroidal  type 
dust  particles  under  the  influence  of  PR  drag,  n(a)  is  quite  simple  (Dohnanyi,  1978) 

.  ^       1 
n[a)  ~  -  .  (4.4) 

a 

For  cometary  type  dust  particles,  because  of  their  high  eccentricities,  the  above 
equation  is  no  longer  valid.  There  is  no  easy  way  that  we  can  find  an  analytic 
expression  for  the  number  density  function.  Therefore  we  have  decided  to  use  the 
empirical  results  from  our  dynamical  study  to  handle  this  problem.  In  our  numerical 
integration,  we  have  sent  out  waves  of  particles  continually  (once  every  300  years) 
from  the  past  until  most  of  the  particles  from  the  first  wave  have  reached  0.1  AU. 
Thus,  the  time-independent  number  density  of  the  particles  between  0. 1  AU  and  the 
source  can  be  obtained  by  analyzing  the  number  of  orbits  per  unit  semimajor  axis 
range  from  the  end  results  of  all  of  our  numerical  integrations.  We  plot  this  result  in 
Figure  4.13.  Apparently  it  is  quite  different  from  the  distribution  of  asteroidal  dust 
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particles,  as  described  by  Equation  4.4.  Because  particles  at  different  semimajor  axis 
ranges  have  different  distributions  in  their  orbital  elements,  it  is  crucial  to  have  the 
correct  expression  for  n(a).  We  have  also  done  a  test  for  our  empirical  method  here. 
We  calculate  the  empirical  number  of  orbits  per  unit  semimajor  axis  for  Eos  dust 
particles  and  compare  it  with  Equation  4.4.  They  agree  very  well  (Figure  4.14). 

The  next  step  is  to  get  the  forced  and  proper  elements  at  different  semimajor 
axis  ranges.  As  an  example,  the  forced  and  proper  eccentricities  from  the  second 
model  are  shown  in  Figures  4.15  and  4.16.  These  quantities  are  well-defined  when 
compared  with  the  results  from  asteroidal  dust  particles.  Thus,  instead  of  using  the 
tabular  form  for  eccentricities  and  longitudes  of  pericenters,  as  we  did  for  asteroidal 
dust  particles,  we  use  a  better  way  of  determining  e  and  w  for  cometary  dust  particles. 
What  we  have  done  is  to  fit  a  polynomial  to  the  points  in  each  diagram  and  generate 
the  elements  at  any  given  semimajor  axis  according  to  those  polynomials. 

If  what  we  have  for  the  eccentricity  is  good  news,  then  the  bad  news  is  that 
we  do  not  have  well-defined  forced  inclination  and  node  for  cometary  dust  particles. 
The  problem  is  in  determining  the  forced  inclinations  (about  1  or  2  degrees)  and 
nodes  from  those  hugely  scattered  points  in  the  inclination  space.  The  least  squares 
fit  method  gives  a  probable  error  of  several  degrees  in  forced  inclination.  We  will 
show  in  the  next  chapter  that  we  are  able  to  match  the  latitudes  of  the  peak  fluxes 
around  the  sky  by  using  the  forced  inclination  and  node  from  asteroidal  dust  particles 
alone.  However,  asteroidal  particles  alone  cannot  account  for  the  observed  shape  of 
the  cloud.  Cometary  dust  particles  may  be  important  in  determining  the  shape  of  the 
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cloud  at  high  ecliptic  latitudes.  In  what  follows,  we  use  the  forced  inclination  and 
node  from  asteroidal  dust  particles  for  the  cometary  dust  particles  (these  are  discussed 
in  Chapter  5).  Since  the  forced  inclinations  from  asteroidal  dust  particles  are  within 
the  probable  errors  of  the  forced  inclinations  from  the  cometary  dust  particles,  we 
feel  that  this  may  be  the  best  way  to  handle  the  problem. 

In  the  next  chapter  we  will  show  how  we  created  zodiacal  cloud  models  based 
on  our  dynamical  study  in  this  chapter  and  previous  chapters  and  will  compare  the 
results  with  the  observations. 
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CHAPTER  5 
STRUCTURE  OF  THE  ZODIACAL  CLOUD 


During  its  all  sky  survey  in  1983,  The  Infrared  Astronomical  Satellite  (IRAS) 
mapped  the  whole  sky  at  four  different  wavelengths:  12,  25,  60,  and  100  ^m.  The 
features  of  the  zodiacal  cloud  are  best  seen  in  the  25  fim  waveband  because  of 
the  weak  galactic  background.  In  this  chapter  we  will  concentrate  on  modeling  the 
zodiacal  cloud  in  the  25  fim  waveband.  All  the  observational  data  shown  hereafter 
are  from  25  ^m  waveband  data. 

The  two  most  important  global  features  of  the  zodiacal  cloud  obtained  from  the 
IRAS  observations  are  the  overall  shape  of  the  cloud  and  the  latitudes  of  the  peak 
fluxes  around  the  sky.  Figure  5.1  shows  the  shape  of  the  zodiacal  cloud  at  90° 
elongation  angle  from  the  observation  sequence  SOP406.  Elongation  angle  is  the 
angle  from  the  Sun  to  the  line  of  sight  of  the  observation.  The  latitude  of  the  peak 
flux  is  close  to  the  ecliptic,  though  not  exactly  in  the  ecliptic.  The  flux  decreases 
towards  high  ecliptic  latitudes.  There  are  several  tiny  bumps  near  the  peak  and  at 
±  10°  ecliptic  latitudes.  These  are  the  so-called  solar  system  dust  bands.  They  are 
caused  by  micron-sized  dust  particles  produced  from  several  prominent  Hirayama 
families  in  the  asteroid  belt  (Dermott  et  al.,  1985,  1993).  The  feature  at  -60°  is 
due  to  the  Galaxy.  Since  we  are  interested  in  interpreting  the  global  structure  of  the 
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zodiacal  cloud,  we  need  to  remove  these  tiny  features  from  the  observations.  The 
filtering  process  is  the  same  as  that  described  in  Dermott  et  al.  (1986). 

Figure  5.2  shows  the  latitudes  of  the  peak  fluxes  at  90°  elongation  angle  from 
the  leading  (open  circles)  and  trailing  (filled  circles)  directions  obtained  by  IRAS 
from  different  positions  along  the  Earth's  orbit  (Dermott  et  al.,  1992).  IRAS  moves 
around  the  Earth  along  the  terminator  with  constant  elongation  angle  during  each 
revolution.  The  leading  direction  is  in  the  direction  of  the  Earth's  orbital  motion, 
while  the  trailing  direction  is  in  the  opposite  sense.  The  two  sine  curves  are  the  best 
fit  results.  The  two  vertical  lines  indicate  where  the  effective  plane  of  symmetry  of 
the  zodiacal  cloud  intersects  the  ecliptic.  Only  at  these  two  nodes  can  the  latitudes 
of  the  peak  fluxes  from  the  leading  and  trailing  directions  have  the  same  magnitude, 
one  above  and  the  other  below  the  ecliptic. 

The  overall  shape  of  the  cloud  and  the  latitudes  of  the  peak  fluxes  are  the  most 
important  observational  constraints  on  the  zodiacal  cloud  model.  They  are  strongly 
related  to  the  dynamical  characteristics  and  evolution  of  the  dust  particles.  Only  that 
model  which  includes  the  correct  particle  dynamics  can  account  for  these  features.  In 
the  following  sections  we  apply  the  dynamical  results  from  the  previous  chapters  to  a 
three  dimensional  model,  SIMUL,  to  construct  a  model  zodiacal  cloud  and  compare 
the  model  with  the  observations.  We  present  the  all-asteroidal  zodiacal  cloud  model 
in  the  first  section.  In  the  second  section  we  show  two  zodiacal  cloud  models  based 
on  cometary  particles  alone  and  one  model  with  combinations  from  both  asteroidal 
and  cometary  particles.  We  give  the  discussion  in  the  third  section. 
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Contribution  from  Asteroidal  Dust  Particles 

In  order  to  construct  a  zodiacal  cloud  model  to  compare  with  the  observational 
data,  we  need  a  model  which  can  generate  the  observable  flux  from  given  distributions 
of  the  orbital  elements  of  the  dust  particles.  The  model  is  called  SIMUL  (Dermott  et 
al.,  1993,  Xu  et  al.,  1993).  What  SIMUL  does  can  be  summarized  as  the  following: 
(1)  generates  millions  of  orbits  according  to  our  dynamical  results,  (2)  defines  a  total 
effective  cross  section  area  and  divides  it  by  the  number  of  orbits,  (3)  spreads  the  area 
along  each  orbit  according  to  Kepler's  law,  (4)  divides  the  three  dimensional  space 
into  millions  of  cells,  (5)  calculates,  for  each  cell,  the  contribution  from  each  orbit, 
and  (6)  produces  a  three  dimensional  cell  table  which  has  the  spatial  distribution 
of  the  effective  cross  section  area.  After  the  table  has  been  generated,  we  integrate 
along  the  line  of  sight  to  obtain  the  effective  flux  in  that  direction.  In  this  way,  we 
are  able  to  mimic  the  viewing  geometry  of  IRAS  and  reproduce  the  observed  flux. 

From  the  results  in  Chapter  3,  we  have  the  forced  elements  as  functions  of 
semimajor  axis  for  9  /xm  diameter  dust  particles.  In  order  to  calculate  the  contribution 
from  all  asteroids  to  the  zodiacal  cloud  we  need  a  bias-free  sample  of  asteroids.  We 
have  chosen  a  sample  of  845  asteroids  between  2.1  AU  and  3.3  AU  with  absolute 
magnitudes  smaller  than  10  as  our  sample.  This  sample  is  bias-free  in  the  sense  that 
it  eliminates  the  observational  bias  against  fainter  asteroids  towards  the  outer  edge 
of  the  belt  (D.  Durda,  personal  communication,  1993).  The  proper  elements  of  these 
asteroids  are  from  Milani  and  Knezevic  (1990).  We  also  assume  that  the  asteroids 
in  this  sample  have  the  same  distribution  of  orbital  elements  as  the  initial  orbits 
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of  the  dust  particles  in  the  cloud,  apart  from  those  changes  produced  by  radiation 

pressure  on  release  from  the  parent  body.  We  then  use  SEMUL  to  calculate  the 
contribution  from  this  bias-free  sample  of  orbital  elements.  The  actual  computation 
in  this  section  was  performed  by  Dr.   Yulin  Xu. 

Comparisons  of  the  shape  of  the  zodiacal  cloud  obtained  from  the  observations 
and  from  the  model  for  two  viewing  geometries  are  shown  in  Figures  5.3  and  5.4.  The 
bold  lines  are  the  IRAS  observations  while  the  thin  lines  represent  the  results  from 
the  models.  The  SOP  number  indicates  a  particular  IRAS  observation.  EA  stands 
for  elongation  angle.  L  and  T  stand  for  leading  direction  and  trailing  direction, 
respectively.  The  comparison  between  the  observations  and  the  results  from  this 
all-asteroidal  model  is  not  at  all  good.  There  is  a  mismatch  between  these  two  sets 
of  data.  At  high  ecliptic  latitude,  the  observations  show  more  flux  than  the  asteroidal 
particles  can  produce.  This  discrepancy  is  seen  at  all  elongation  angles  and  at  all 
longitudes  where  the  observations  were  made. 

The  latitudes  of  the  peak  fluxes  around  the  sky  from  the  all-asteroid  model, 
however,  fit  the  observations  very  well.  Figure  5.5  shows  the  comparison  for  90° 
elongation  angle  around  the  sky.  Thin  lines  are  the  results  from  the  all-asteroid 
model.  Bold  lines  are  the  best  fit  curves  from  the  IRAS  observations.  The  differences 
at  the  ascending  node  are  0°.08  and  12°  .7  for  the  inclination  and  longitude  of  the 
ascending  node,  respectively.  The  differences  at  the  descending  node  are  0°.01  and 
6°  .5  for  the  inclination  and  longitude  of  the  descending  node,  respectively. 
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Contribution  from  Encke-type  Dust  Particles 

In  Chapter  4  we  have  two  models  for  Encke-type  dust  particles  with  different 
initial  inclination  distributions.  The  results  from  the  first  model  are  presented  first. 
SIMUL,  again,  is  used  to  calculate  the  flux  from  the  dust  particles.  In  Figure  5.6 
we  show  the  shape  of  the  zodiacal  cloud  at  SOP392,  97°  elongation  angle,  leading 
direction,  from  the  model  (thin  curves)  and  from  the  observation  (bold  curves). 
There  are  three  panels  in  the  diagram,  the  top  one  being  the  contribution  from  the 
all-asteroid  model.  The  second  one  is  the  contribution  from  the  first  Encke-type 
dust  particle  model  while  the  third  panel  is  a  combination  of  27%  from  the  all- 
asteroid  model  and  73%  from  the  first  Encke-type  dust  model.  It  is  obvious  that 
the  observation  cannot  be  matched  unless  we  use  the  right  combination  of  asteroidal 
and  cometary  particles.  Even  though  the  match  in  the  third  panel  is  not  perfect,  it 
is  a  good  match. 

In  Figure  5.7  we  show  a  similar  diagram  for  the  shape  of  the  zodiacal  cloud  at 
SOP406,  90°  elongation  angle,  leading  direction.  Again,  the  best  combination  seems 
to  be  73%  cometary  particles  and  27%  asteroidal  particles.  We  have  found  that  this 
combination  is  independent  of  the  position  of  observation  and  the  elongation  angle. 

The  latitudes  of  the  peak  fluxes  around  the  sky  from  the  combination  of  27% 
asteroid  particles  and  73%  cometary  particles  (from  the  first  model)  are  shown  in 
Figure  5.8.  These  are  the  data  at  90°  elongation  angle  around  the  sky.  When  we 
compare  this  with  the  observations,  the  differences  at  the  ascending  node  are  0°.37 
and  27°. 4  for  the  inclination  and  longitude  of  the  ascending  node,  respectively. 
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The  differences  at  the  descending  node  are  0°.28  and  21°. 6  for  the  inclination  and 
longitude  of  the  descending  node,  respectively.  Clearly,  this  is  not  as  good  as  the 
all-asteroidal  particle  model. 

The  zodiacal  cloud  constructed  from  the  second  cometary  particle  model  is  not 
as  good  as  that  using  the  first  one.  The  comparison  between  the  observations  and  the 
all-comet  model  at  two  different  SOP  positions  are  shown  in  Figure  5.9.  North  of  the 
ecliptic,  there  is  more  flux  from  the  model  than  the  observations.  However,  south  of 
the  ecliptic,  the  model  cannot  even  produce  enough  flux  to  match  the  observations. 
Thus,  it  is  impossible  to  match  the  shape  of  the  zodiacal  cloud  by  combining  the 
asteroid  model  and  this  cometary  particle  model.  The  difference  between  these  two 
cometary  models  will  be  discussed  in  the  next  section. 

Discussion 

Based  on  our  dynamical  study  and  modeling  results,  we  can  draw  several 
conclusions  about  the  structure  of  the  zodiacal  cloud: 

1.  Single-sized  asteroidal  dust  particles  alone  match  very  well  the  latitudes  of  the 
peak  fluxes  around  the  sky. 

2.  Single- sized  asteroidal  dust  particles  alone  cannot  come  up  with  enough  flux  at 
high  ecliptic  latitudes. 

3.  Cometary  dust  particles  cannot  form  a  dust  band. 

4.  Cometary  dust  particles  can  provide  more  flux  at  high  ecliptic  latitudes. 


112 
5.    By  choosing  a  suitable  initial  inclination  for  cometary  dust  particles,  we  have 
found  a  good  fit  for  the  shape  of  zodiacal  cloud,  with  a  73%  contribution  coming 
from  cometary  particles  and  the  remaining  27%  coming  from  asteroidal  particles, 
but  the  latitudes  of  the  peak  fluxes  around  the  sky  are  not  well  reproduced. 

In  coming  up  with  these  conclusions,  we  have  been  using  a  single-sized  particle 
model  and  we  have  not  considered  collisions  between  particles  while  they  are 
spiralling  towards  the  Sun.  Most  of  the  9  //m  particles  produced  at  the  asteroid 
belt  have  collisional  life  time  much  longer  than  the  PR  drag  life  time.  More  than 
95%  of  them  make  it  to  1  AU  without  being  destroyed  by  a  collision,  while  50  % 
of  50  fim  diameter  particles  will  be  coUisionally  shattered  before  they  reach  lAU 
(Gustafson,  1994). 

Conclusion  1  certainly  is  a  significant  accomplishment  for  our  dynamical  study 
and  modeling  process.  Conclusion  2  has  posed  a  challenge  to  the  all-asteroid  model. 
Is  it  possible  that  by  including  the  collisional  process  of  dust  particles  we  may  find 
a  solution  to  this  problem?  This  is  a  tough  question  which  cannot  be  answered 
easily.  Some  arguments  can  be  addressed  concerning  this  alternative.  First  of  all, 
to  include  the  collisional  process  of  the  particles  is  equivalent  to  adding  additional 
source  terms  between  the  asteroid  belt  and  Sun.  So  instead  of  having  the  number 
density,  n{r)  ~  r~\  we  have  n{r)  ~  r~'",  where  m  is  larger  than  one.  Satellite 
experimental  data  (see  the  review  papers  by  Dohnanyi,  1978,  Leinert  and  Griin, 
1990)  indicate  m  ~  1.3  with  some  uncertainties.  With  the  model  (SIMUL),  the 
effect  of  additional  sources  on  the  shape  of  the  flux  at  high  ecliptic  latitudes  can  be 
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tested  easily  by  changing  the  index  from  1  to  ,  say,  2  and  seeing  if  the  flux  at  high 
ecliptic  latitude  can  be  brought  up  to  match  the  observations.  The  answer  is  no.  This 
is  understandable.  At  lower  ecliptic  latitude  we  will  gain  the  additional  flux  all  the 
way  to  the  asteroidal  belt.  On  the  other  hand,  at  high  ecliptic  latitude,  we  only  gain 
the  additional  flux  up  to  a  certain  distance  from  the  Earth.  Thus,  relatively  speaking, 
we  gain  more  flux  at  low  ecliptic  latitude  than  at  high  ecliptic  latitude.  Unless  the 
index  is  so  large  that  there  is  a  sharp  increase  near  the  Earth  (this  is  not  impossible), 
the  peak  will  only  become  sharper. 

Another  argument  concerning  collisional  products  from  large  particles  is  the 
change  in  their  forced  elements.  Large  particles  (e.g.  50  ^m  in  diameter)  have  a 
very  slow  drag  rate  and  will  pick  up  large  forced  inclinations  on  their  way  to  the 
Sun  when  they  pass  through  secular  resonance  regions.  This  means  when  small  dust 
particles  are  produced  from  these  large  particles,  somewhere  between  the  Earth  and 
the  asteroid  belt,  they  will  start  with  some  very  different  forced  elements  than  their 
parent  bodies'  at  the  asteroid  belt.  The  evolution  and  distribution  of  these  particles 
will  be  quite  different  from  those  with  the  same  size  but  produced  from  the  asteroid 
belt.  When  we  add  the  contribution  from  these  particles  to  our  model,  will  their 
forced  elements  change  the  plane  of  symmetry,  i.e.,  the  latitudes  of  the  peak  fluxes 
around  the  sky?  Will  we  still  be  able  to  obtain  a  good  fit  between  the  model  and 
the  observations?  This  question  and  those  we  have  just  mentioned  in  the  previous 
paragraph  are  not  easy  questions  to  answer.  Until  we  include  the  collision  process 
in  the  dynamical  model,  we  cannot  know  what  the  result  is. 
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Another  way  to  increase  flux  at  high  ecliptic  latitude  for  the  all-asteroid  model 
is  to  add  more  high  inclination  asteroids  to  the  sample  we  used  to  generate  the 
zodiacal  background  in  the  previous  section.  This  is  possible.  There  may  exist  an 
observational  bias  against  high  inclination  asteroids.  In  other  words,  the  sample  we 
use  to  represent  the  whole  asteroid  population  may  not  be  bias  free  in  its  proper 
inclination  distribution.  But  it  is  unlikely  that  the  bias  against  high  inclination 
asteroids  in  the  observations  is  so  large  that  when  we  put  the  real  bias  free  sample 
in  our  calculation  we  will  be  able  to  obtain  enough  flux  at  high  ecliptic  latitude. 

Comets,  due  to  their  high  eccentricities,  certainly  provide  a  possible  way  to  add 
flux  to  the  high  ecliptic  latitude.  However,  from  our  two  models  we  can  see  that 
only  with  the  right  distribution  in  initial  inclinations  can  we  come  up  with  enough 
flux  at  high  ecliptic  latitude.  The  distribution  in  inclinations  for  all  Jupiter  Family 
(JF)  short  period  comets  is  shown  in  Figure  5.10.  These  are  the  latest  available 
observational  data  (Marsden  and  Williams,  1992).  We  don't  know  whether  this 
represents  a  time-independent  population  or  not.  There  may  also  exist  a  strong  bias 
against  high  inclination  comets.  Even  if  this  is  an  equilibrium  population,  we  still 
don't  know  how  many  of  them  actually  contribute  to  the  zodiacal  cloud  and  what  are 
their  relative  contributions.  What  we  do  know  is  that  Encke  itself  is  very  unlikely 
to  be  the  major  source  which  supplies  dust  to  the  zodiacal  cloud  because  of  the 
result  from  our  second  cometary  model.  If  there  are  more  high  inclination  JF  comets 
depositing  dust  particles  into  the  interplanetary  space,  then  we  may  have  less  than 
70%  of  the  zodiacal  cloud  coming  from  comets.  If  there  are  too  many  low  inclination 
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JF  comets  that  deposit  dust  to  the  zodiacal  cloud,  then  we  cannot  even  use  cometary 
dust  particles  to  bring  up  enough  flux  at  high  ecliptic  latitude. 

Dust  particles  produced  from  a  comet  with  its  apocenter  outside  the  orbit  of 
Jupiter  will  have  their  orbits  cross  that  of  Jupiter.  Most  of  the  particles  will  be 
gravitationally  scattered  by  Jupiter.  Thus,  comets  with  their  apocenters  inside  the 
orbit  of  Jupiter  are  likely  to  be  the  most  effective  sources  to  deposit  dust  particles 
into  the  zodiacal  cloud.  In  Figure  5.11  we  show  the  distribution  in  eccentricity  and 
inclination  for  those  JF  comets  which  have  apocenters  inside  the  orbit  of  Jupiter  and 
eccentricities  larger  than  0.5  (so  that  the  inclinations  of  the  dust  particles  produced 
from  their  components  can  go  through  large  enough  variations).  There  are  12  of  them. 
Eight  of  them  have  inclinations  larger  than  9°  while  only  4  of  them  have  inclinations 
smaller  than  6°.  Of  course  this  may  not  be  a  time  independent  population.  But  if 
it  is,  and  if  they  do  contribute  equally  to  the  zodiacal  cloud,  we  would  expect  the 
dust  produced  from  this  population  to  have  an  inclination  distribution  concentrated 
towards  high  inclination  similar  to  that  in  Figure  4.11.  If  such  is  the  case,  then  we 
can  form  a  zodiacal  cloud  with  the  major  component  coming  from  these  comets. 

The  latitudes  of  the  peak  fluxes  around  the  sky  are  not  well  reproduced  when 
we  combine  73%  of  the  first  cometary  particles  with  27%  of  the  asteroid  particles  to 
form  the  zodiacal  cloud.  This  is,  again,  related  to  the  actual  source  of  the  cometary 
particles.  It  is  true  that  we  apply  the  same  forced  inclination  and  node  from  the 
asteroidal  model  to  the  cometary  model.  However,  what  we  see  from  the  observations 
is  an  effective  plane  of  symmetry.  It  is  a  result  of  integrating  along  the  line  of  sight 
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through  the  warped  plane  of  symmetry  of  the  cloud.  Because  cometary  particles  have 
different  distributions  in  their  eccentricities  and  inclinations  than  those  of  asteroid 
particles,  it  is  not  surprising  that  the  combined  result  deviates  from  the  all-asteroid 
model.  If  we  are  able  to  determine  the  real  source  of  cometary  particles,  then,  by 
studying  statistically  enough  dust  particles  from  the  source,  we  may  be  able  to  find 
the  right  forced  inclination  and  node  for  cometary  particles  and  match  the  latitudes 
of  the  peak  fluxes  as  well. 

Near  Earth  asteroids  (NEA)  are  another  possible  source  which  could  contribute 
dust  particles  to  the  zodiacal  cloud.  They  produce  dust  particles  through  collisions 
between  themselves  and  with  main  belt  asteroids.  Due  to  their  similar  distributions  in 
eccentricities  and  inclinations  with  short  period  comets,  dynamically  speaking  NEA 
and  Jupiter  Family  comets  are  similar  groups,  if  not  identical.  NEA  are  capable  of 
producing  the  same  kind  of  dust  that  JF  comets  produce.  The  bad  news  is,  all  the 
problems  that  JF  comets  have  apply  to  NEA,  also.  We  do  not  know  whether  the 
observed  NEA  represents  a  bias-free  population  or  not.  Nor  do  we  know  their  relative 
contribution.  It  may  also  be  possible  that  both  JF  comets  and  NEA  contribute  to 
the  zodiacal  cloud. 

Another  possible  way  to  account  for  the  flux  at  high  ecliptic  latitude  is  to  add  a 
spheroidal-shaped  flux  to  the  result  from  the  asteroid  model.  From  the  mathematical 
point  of  view,  this  is  a  model  which  can  fit  the  observations.  However,  the  source 
of  this  spheroidal-shaped  flux  is  unknown.  Maybe  it  is  due  to  the  contribution  from 
interstellar  grains.  But  there  are  too  many  unknowns  to  this  alternative. 
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From  our  dynamical  study  in  this  dissertation,  we  have  obtained  an  understanding 
of  the  orbital  evolution  of  both  asteroidal  and  cometary  dust  particles.  We  have 
examined  all  the  effects  that  they  can  face  in  the  real  solar  system.  By  using  our 
dynamical  results  from  particles  of  one  single  size  to  construct  several  zodiacal  cloud 
models,  we  have  concluded  that  a  zodiacal  cloud  constructed  from  different  kinds  of 
particles  (or  with  different  contributions  from  different  kinds  of  particles)  will  have 
quite  different  characteristics.  We  show  that,  by  the  right  combination,  it  is  possible 
to  build  a  model  which  can  fit  the  overall  shape  of  the  zodiacal  cloud  from  the  IRAS 
observations.  This  may  be  a  major  step  towards  revealing  the  true  structure  of  the 
zodiacal  cloud.  The  next  step  will  be  to  include  particles  with  a  range  of  sizes  in  the 
zodiacal  cloud  model  and  to  include  collisional  effects  in  the  dynamical  evolution 
of  dust  particles.  In  addition,  the  fundamental  question  needed  to  be  answered  next 
is,  what  are  the  exact  sources  of  the  dust  particles?  There  are  numerous  as-yet 
unanswered  questions  about  the  main  belt  asteroids,  JF  comets,  and  NEA:  is  there 
an  observational  bias  against  high  inclination  asteroids?  how  large  is  this  bias?  are 
JF  comets  a  stable  population?  how  many  of  them  contribute  to  the  zodiacal  cloud? 
do  NEA  produce  a  significant  amount  of  dust  particles?  Until  these  questions  are 
answered,  the  actual  sources  of  the  zodiacal  cloud  cannot  be  determined.  Nor  can 
we  determine  the  real  constituents  of  the  zodiacal  cloud. 
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